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Abstract. This review is concerned with the gravitational self-force acting on 
a mass particle in orbit around a large black hole. Renewed interest in this 
old problem is driven by the prospects of detecting gravitational waves from 
strongly gravitating binaries with extreme mass ratios. We begin here with a 
summary of recent advances in the theory of gravitational self-interaction in 
curved spacetime, and proceed to survey some of the ideas and computational 
strategies devised for implementing this theory in the case of a particle orbiting 
a Kerr black hole. We review in detail two of these methods: (i) the standard 
mode-sum method, in which the metric perturbation is regularized mode-by-mode 
in a multipole decomposition, and (ii) m-mode regularization, whereby individual 
azimuthal modes of the metric perturbation are regularized in 2+1 dimensions. 
We discuss several practical issues that arise, including the choice of gauge, 
the numerical representation of the particle singularity, and how high-frequency 
contributions near the particle are dealt with in frequency-domain calculations. 
As an example of a full end-to-end implementation of the mode-sum method, we 
discuss the computation of the gravitational self-force for eccentric geodesic orbits 
in Schwarzschild, using a direct integration of the Lorenz-gauge perturbation 
equations in the time domain. With the computational framework now in place, 
researchers have recently turned to explore the physical consequences of the 
gravitational self force; we will describe some preliminary results in this area. An 
appendix to this review presents, for the first time, a detailed derivation of the 
"regularization parameters" necessary for implementing the mode-sum method in 
Kerr spacetime. 
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1. Introduction 



Context and motivation. Within Newtonian theory, the gravitational two-body 
problem, in its basic form, is readily solvable. An isolated system of two 
gravitationally-bound point masses admits two conserved integrals — the energy 
and angular momentum — and the resulting motion is precisely periodic. The 
corresponding general- relativistic problem is radically — and notoriously — more 
difficult. In General Relativity (GR), the orbits in a bound binary are never periodic: 
Gravitational radiation continually removes energy and angular momentum from the 
system, and the radiation back-reaction gradually drives the two objects tighter 
together until they eventually merge. Furthermore, in GR one cannot usually work 
consistently with pointlike mass particles (as we discuss later), and the simplest and 
most universal two-body problem becomes that of a binary black hole system. 
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In general, the description of the nonhnear radiative dynamics of a black hole 
binary entails a full Numerical-Relativistic (NR) treatment. There has been a much 
celebrated advance in NR research in the past few years [1381 H51 H] , with NR codes 
now capable of tracking the complicated nonlinear evolution of a binary black hole 
spacetime during the final stages of the merger. However, NR methods become less 
efficient when the two black holes are far apart, or when one of the components is 
much heavier than the other. Each of these two regimes features two greatly different 
lengthscales, which is not easily accommodated in a NR framework [76]. Fortunately, 
the occurrence of two separate lengthscales also restores a sense of "pointlikeness" 
(at least approximately), which allows for simpler, perturbative treatments. In 
the first regime — at sufficiently large separations — the dynamics is best analyzed 
using post-Newtonian (PN) methods (3^, wherein GR corrections to the Newtonian 
dynamics (accounting for radiation reaction, the objects' internal structure, etc.) are 
incorporated into the equations of motion order by order in the binary separation. 

The second, so-called extreme mass-ratio regime — which will concern us in this 
review — is most naturally explored within the framework of black hole perturbation 
theory. Here the "zeroth-order" configuration is that of a test particle (the lighter black 
hole) moving along a geodesic of the fixed background spacetime of the large black 
hole. This can then serve as a basis for a systematic expansion, wherein corrections due 
to the finite mass of the small object (and due possibly also to its internal structure) 
are included order by order in the small mass ratio. At first order in the mass ratio, the 
gravitational field of the small object is a linear perturbation of the background black 
hole geometry. The back reaction from this perturbation gives rise to a gravitational 
self force (SF) which gradually diverts the small object from its geodesic motion. In 
this picture, it is the SF that is responsible (in particular) for the radiative decay 
of the orbit. In principle, knowledge of the SF (along with the metric perturbation) 
forms a complete picture of the orbital dynamics at linear order in the small mass 
ratio. One therefore hopes that a calculation of the SF would facilitate a faithful, if 
only approximate description of the orbital dynamics in binary systems with small 
mass ratios. 

The determination of the SF in curved spacetime is an old problem in 
mathematical relativity, made very relevant following recent developments in 
gravitational-wave research. Interest in this problem was renewed in the mid 1990s, 
when it was first proposed that the planned space-based gravitational wave detector 
LISA (the Laser Interferometer Space Antenna [2]) could observe signals from the 
inspiral of compact objects (white dwarfs, neutron stars, or stellar-mass black holes) 
into massive black holes in galactic nuclei. It is now believed that LISA should 
be able to detect tens to thousands such events [70], out to cosmological distances 
{z ~ 1) [68j — depending on the astrophysical rates (inspirals per galaxy per year), 
which are still highly uncertain^ Dubbed EMRIs (Extreme Mass Ratio Inspirals), 
these potential sources became key targets for LISA due to their unique facility as 
precision probes of strong-field gravity. A typical LISA EMRI will spend the last few 
years of inspiral in a very tight orbit around the massive hole, emitting some 10^-10^ 
gravitational wavecycles wholly within the LISA frequency band. These complicated 
waveforms carry extremely accurate information about the physical parameters of 
the inspiral system, as well as a detailed map of the spacetime geometry around the 

X In fact, if inspiral rates are near the higher end of their estimated range, the unresolvable stochastic 
background from thousands weak inspirals could dominate the noise budget of LISA at its most 
sensitive frequency band. |15| 
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massive hole. The potential scientific implications are broad and far-reaching, ranging 
from astrophysics to cosmology and to fundamental theory [m[5ni[TF[ [51 [7TllT^ll48j . 

The odds are that even the strongest detected EMRI signals will be buried deeply 
inside LISA's instrumental noise. However, since an EMRI signal is detectable over a 
long time, it can be extracted with significant signal-to-noise ratio (of order ~ 100 for 
the strongest sources) using matched filtering [14j . This, however, requires precise 
theoretical templates of EMRI waveforms across the relevant parameter space of 
inspirals — indeed, the detection statistics cited above assumes that such templates 
were at hand by the time LISA flies. It is likely that other, non-template-based 
data analysis techniques could be used to detect some of the brightest EMRIs 
[72l[8]. Nonetheless, it remains true that an accurate parameter extraction, crucial for 
exploiting the full scientific value of the EMRI signal, will rely on (or be restricted by) 
the availability of accurate and faithful theoretical templates of the inspiral waveforms. 

The theoretical challenge is derived from the astrophysical specifications of the 
LISA-relevant EMRIs. With its peak sensitivity at a few mHz, LISA will observe 
inspirals into galactic (likely Kerr) black holes with masses in the range ~ 5 x 10^- 
5 X 10^ Mq; the inspiraling objects (compact stars or stellar-mass black holes) may 
have masses in the range O.5-5OM0, giving mass ratios of 10~^-10~® — well within the 
"extreme mass ratio" domain of the two-body problem. Inspirals are not expected 
to have any preferred orientation with respect to the central hole's spin direction, 
and orbits may remain quite eccentric throughout the inspiral [M]. It is not likely 
that interaction with a possible accretion disk around the massive hole would play 
a dominant role in the orbital dynamics [3T]. In a typical LISA-band EMRI — a 
IOMq/IO^Mq system — gravitational radiation reaction drives the orbital decay over 
a timescale of months. More importantly, it affects the phasing of the inspiral orbit 
over mere hours. It is therefore clear that a useful model of the long-term orbital 
phase evolution in a LISA-relevant system ought to incorporate properly the effect of 
radiation reaction. 

The above translates, at first approximation, to a very clean problem in black hole 
perturbation theory: A point mass is set in a generic (eccentric, inclined) strong-field 
orbit around a Kerr black hole of a much larger mass, and one wishes to calculate the 
gravitational waveforms emitted as radiation reaction drives the gradual inspiral up 
until the eventual plunge through the event horizon. While it remains important to 
quantify the effect of higher-order corrections to this picture (e.g., due to the spin of 
the inspiralling object), one expects that, thanks to the extreme mass ratio, the above 
simple setup should provide a good model for astrophysical inspirals. Indeed, from the 
perspective of GR theorists working in the field, much of the appeal of the SF problem 
comes from its unusual dual nature as both an elementary theoretical problem in GR, 
and an exciting problem in contemporary astrophysics. 

Scope and relation to other reviews. Our main aim here is to review the challenges and 
main developments in the program to calculate the gravitational SF in Kerr spacetime. 
Our discussion will be focused mostly on work concerned with the evaluation of the 
SF along a specified, non-evolving orbit (normally taken to be a geodesic of the 
background spacetime); we will not consider here the important question of how orbits 
evolve under the effect of the SF. The strategic approach envisaged here is one in which 
the complete analysis of the orbital evolution is carried out in two separate, consequent 
steps. In the first, preparatory step, one calculates the SF across the entire relevant 
phase space (i.e., obtain the value of the SF as a function of location and velocity. 
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perhaps through interpolation of numerical results). In the second step one then 
uses the SF information to calculate the inspiral orbits of particles with given initial 
conditions. Here we will be concerned only with the first, most crucial step. 

In parallel to the work on SFs described in this review, there is already a 
substantial research effort aimed to formulate a reliable scheme for the orbital 
evolution, assuming the SF has been calculated [1371 l87] . This work incorporates 
techniques from multiple-scale perturbation theory. Another parallel effort aims to 
obtain an approximate model of the orbital evolution and emitted waveforms without 
resorting to the load SF [Ml [Ml El [HSl El [113]. This work is largely based on a 
strategy proposed by Mino |110[ lllli 11451 1112] , in which a time-average measure of 
the rate of change of the geodesic "constants of motion" is calculated from a certain 
"radiative" solution of the perturbation equations (to be described later in this review), 
which accounts for the long-term radiative aspects of the dynamics but neglects some 
of the conservative effects. This "inspiral without SF" approach is reviewed by Mino 
in Ref. [lllj , and also by Tanaka in Ref. [154j . It is well possible that this method will 
prove sufficiently accurate for LISA applications. However, ultimately, its performance 
and accuracy could only be assessed against actual calculations of the full SF. 

The fundamental formulation of the SF in curved spacetime is described in a 
comprehensive Living Review article by Poisson [131j . This is an excellent treatise 
which is both self-contained and pedagogical, and it makes an essential reading for 
anyone wishing to introduce oneself to the field. (The less endurant reader would 
find a slightly abridged version of this review in Ref. [132] : there is also a concise 
introduction in |133j .) It would not be useful to review here at any great detail the 
basic SF theory already covered in |131j . Instead, we merely give a succinct summary 
of the essential formal results, and move on to describe, in the rest of our review, how 
the general formalism is implemented in actual calculations of the gravitational SF in 
Kerr. In this respect, our review starts where Poisson's review ends. We will, however, 
briefly survey a few recent developments in SF theory not covered in Ref. |131j (last 
updated in 2004). 

A good 2005 snapshot of the activity surrounding SF calculations is offered by a 
collection of review articles published in a special issue of Classical and Quantum 
Gravity [IDT. Drasco's review from 2006 [62] explores the utility of the various 
computation methods within the context of the LISA EMRI problem. The website of 
the 12th Capra Meeting on Radiation Reaction in Relativity (Bloomington IN, 2009) 
[1] is an excellent resource of information on current research in the field, with links 
to all talks given in the meeting. The website also includes links to the websites of 
previous meetings in the series (1998-2008). There are several good reviews of EMRI 
astrophysics and the science potential of EMRI detections, including ones by Hopman 
[89] . Amaro-Seoane et al. [3], and Miller et al. jl09j . Finally, an abridged version of 
some of the material in our current review is to appear in Ref. [12] . 

Structure. The overall structure of our presentation will follow the logical route 
of development in SF research: From basic theory (Sec. ^ through to practical 
calculation schemes (Sees. [3]-[4l) and to numerical implementation in Schwarzschild 
(Sees. [IHl]) and Kerr (Sec. [T]), concluding with a discussion of some physical 
consequences (Sec. [8]). 

We begin in Sec. [2] with a brief introduction to the general theory of the 
gravitational SF in curved spacetime, highlighting from the outset the role of gauge 
dependence in this theory. We then summarize (rather than reproduce) the main 
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theoretical developments that over the past twelve years helped establish a robust 
formal framework for SF calculations. Section [3] is a brief survey of the main 
strategies that have been proposed for implementing this formal framework in actual 
calculations, particularly for orbits in Kerr (or Schwarzschild) spacetimes. This section 
also provides a quick-reference catalogue (Tables [THS]) of actual computations of the 
SF carried out so far. Section |4] is a self-contained introduction to the mode-sum 
method, one of the leading techniques for SF calculations. The basic idea is presented 
through an elementary example, followed by a formulation of the method as applied to 
generic orbits in Kerr. In an accompanying appendix we provide, for the first time, a 
full derivation of the regularization parameters necessary for implementing the mode 
sum scheme in Kerr. (The values of these parameters where published in the past [24] 
without a detailed derivation.) 

An essential preliminary step in almost all calculations of the SF involves the 
numerical integration of the relevant perturbation equations sourced by the point 
particle. In Sec. [5] we discuss the practicalities of such calculations, and review some 
of the numerical strategies that were proposed for dealing with the particle singularity 
in both the frequency and time domains. Section [6] then focuses on a particular 
implementation strategy, based on a direct time-domain integration of the Lorenz- 
gauge metric perturbation equations. This approach recently led to a first computation 
of the gravitational SF for eccentric orbits in Schwarzschild, and we show some results 
from this calculation. In Sec.[7|we discuss higher-dimensional alternatives to standard 
mode-sum, focusing on the recently proposed m-mode regularization, which may offer 
a more efficient treatment in the Kerr case. Section |8] reviews initial work aimed to 
understand and quantify the physical, gauge-invariant effects of the gravitational SF. 
Finally, in Sec.lHlwe reflect on recent advances and comment on future directions. 

Notation. We follow here the notation conventions of Misner, Thorne and Wheeler 
[117] . Hence, the metric signature is ( — I— the connection coefficients and 
Riemann tensor are F^^ = \g^"{grj^,^v + gav.ti. - gt^i^.a) and R^x^,^ = F^^^ - F^^_^ -f 
^afjXxu ~ ^a-i^^Xfi^ the Ricci tensor and scalar are Rap = R^afip and R = Ra", and 
the Einstein equations read Ga/s = Rap — \gapR = SttTq,^. We use Greek letters for 
spacetime indices, and adopt standard geometrized units (with c~G~ 1) throughout. 

2. Essential theory 

2.1. Gravitational forces 

The concept of a gravitational force is, of course, fundamentally strange to GR: In 
a purely gravitational system one expects no acceleration and hence no "forces" in 
the ordinary sense. Since the gravitational SF is an example of a force of a purely 
gravitational origin, we begin by explaining what is generally meant by "gravitational 
forces" in our context. 

Consider a smooth region of spacetime with metric gap (which may be, for 
example, the stationary vacuum exterior of a Kerr black hole), and a smooth weak 
gravitational perturbation hap of that spacetime (which may represent, for example, 
an incident gravitational wave). Now consider a test particle of mass /z which is moving 
freely in the perturbed spacetime. Neglecting SF effects, the particle's trajectory will 
be a geodesic of the perturbed spacetime, described, in a given coordinate system a;". 
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by 

dT'^ ^ ''"dT' dT' ' ^ ' 

where r' is an afRne parameter along the trajectory and T'^^^ are the connection 
coefficients associated with the perturbed metric g + h. 

In some occasions, however, it is practicahy useful to reinterpret the particle's 
motion in terms of a trajectory in the background spacetime g. Under this 
interpretation, the trajectory (in g) is no longer geodesic; rather, the particle 
experiences "an external gravitational force" , which is exerted by the perturbation 
h. This (fictitious) force is defined through Newton's second law as 

_ 

^grav-M(^ dr J' ^ ' 

where t is an affine parameter in the background metric g, and F^^ are the connection 
coefficients associated with g. Note that, in this non-covariant description, h and 
r are treated as tensor fields in g, and similarly the force F^^v four-velocity 
u" = dx'^/dr are defined as vectors in g. The tensorial indices of all these quantities 
are raised and lowered using the background metric g. 

To express i^^av terms of the linear perturbation field hafs, we use d/dr = 
{dr' /dT)d/dT' in Eq. Q and introduce AT^^^ = F^"^ - F;^^. By virtue of Eq. ^ this 
gives 

Fg^.,, = -/.AF>V + C«", (3) 
with C, = ii{dT/dT'){d'^T' /dr^). From its definition in Eq. the force F^^v must 
be perpendicular to the four velocity w". Hence, projecting i^^^v orthogonally to u" 
keeps it unchanged, and we may use this fact to dispose of the term Cu" in Eq. 

^g"av = "M(^A + w"ua)AF;:,U^U^ (4) 

Finally, expressing AF in terms of h (keeping only terms that are linear in h) we 
obtain 

Fg",a, = -l^g"^ + u"u^) {V,hx^ + V^/iA. - Va V) = mV^'^^/i/j^, (5) 

where Vq denotes covariant differentiation with respect to the background metric g. 
The differential operator V"''''' determines the gravitational force exerted by any given 
external perturbation; it is given explicitly by 

\/°'P'r = i (c,"*m/3 - 2g"'3y5 _ ^0^/3^5^ ^Ty^. (6) 
Later we will often work with the trace-reversed metric perturbation, 

haf3 = hap — -gapg^^h^^. (7) 

In terms of h^p, Eq. ([5]) becomes 

F«.a, = iiV'^^-'hp^, (8) 

with 

ya/37 ^ 1 [2g"^u'^u''' - ig'^^u'^u^ - 2u°'u^u'u^ + u^g^^u^ + g"^g'^^) Vs- (9) 
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Formally, the above interpretation of the motion requires a suitable procedure for 
mapping the physical trajectory from the full spacetime g + h onto the background 
spacetime g. In the above discussion we adopted the following procedure: First, we 
assume that g -\- h and g were covered with two coordinate meshes that are "similar" 
in the sense that, at the limit /i ^ 0, any given physical event would attain the same 
coordinate values in both spacetimes. Then, we identify each event a;" u\ g + h with 
an event having the same coordinate value in g. This, in particular, results in a 
projection of the trajectory in g + /i onto one in g. 

Of course, there is not just one way of specifying the coordinate systems in the 
two spacetimes, which leads to an ambiguity in the projected trajectory: Different 
coordinate choices would, in general, give rise to different projected trajectories in g, 
showing different accelerations and hence interpreted as being under the influence of 
different gravitational forces. This, of course, is nothing but an example of the usual 
gauge ambiguity intrinsic to perturbation theory in GR. The gravitational force, just 
like the metric perturbation itself, is gauge dependent. 

It is straightforward to write down a gauge transformation law for -F^^av ■ Consider 
a small gauge displacement 

x^'^x^'-e, (10) 

where the magnitude of ^'^ is assumed to scale like that of the external perturbation 
hai3- Under the perturbation transforms as hap hap + 5(^^)hap, where 

<^(5)^a/3 = VaC/3 + V/3$a. (11) 

From Eq. ([5|), this will induce a change in the gravitational force, given by 

(Terms arising from the gauge transformation of V"''''' are quadratic in the magnitude 
of the perturbation and we neglect them here.) Substituting from Eq. ()11|) in Eq. 
(|12p and using the commutation relation V^Vj/^a — V,yV^^a = £,\R^aiyij.i one readily 
arrives at 

-^(O^g^rav = -M + W'u^) L + R''^.x.u^^eu''] , (13) 

where overdots denote covariant derivatives with respect to the afhne parameter r 
along the background trajectory. This is the general gauge transformation formula for 
external gravitation forces 



2.2. Gravitational self force 

To devise a theory of the gravitational SF, one might be tempted to simply interpret 
it as an example of a gravitational force of the type discussed above, with the source 
of the metric perturbation now being the particle itself. This naive interpretation 
would be problematic, for several reasons. First, the physical perturbation due to the 
particle (a retarded solution of the linearized Einstein equations, h'^^p) is singular at 
the location of the particle, and the statement that the particle follows a geodesic of 

§ One might notice that Eq. JTSj reproduces the geodesic deviation equation if one sets the left- 
hand side to zero and reinterprets (^'^ as the displacement vector connecting two adjacent geodesies 
in g. Indeed, a vanishing 6(^qF^^^^ implies that the two adjacent trajectories (the original projected 
trajectory and its gauge-transformed counterpart) have the same acceleration in g, in which case the 
displacement vector between them is known to satisfy the geodesic deviation equation. 
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g + is therefore physically meaningless. Obviously, trying to apply Eq. ^ [or ([5])] 
with the external perturbation replaced with the self-perturbation h"^^ would yield a 
singular, and hence meaningless result. 

Second (and relatedly), since we are now considering the self-gravity of the 
particle (it is no longer a test particle), we must make a mathematical sense of 
its being "pointlike" . This is not a trivial matter to address in curved spacetime. 
Mathematically, the usual delta- function representation of a point particle stress- 
energy is known to be inconsistent with the nonlinearity of the full Einstein equations 
[7^ I161j . A familiar physical manifestation of this is that one cannot squeeze a 
finite amount of mass to a point without creating a black hole (of a finite size) . The 
mathematical consistency of a delta-function source is restored in the linear theory, 
but it remains a challenging task to understand how the notion of a point particle 
might emerge (rather than be pre-assumed) from a suitable limiting procedure. 

Third, the SF is conceptually different from the external forces discussed above, 
in that the latter are, in truth, just fictitious forces resulting from our insistence to 
artificially split the physical spacetime into a background and a perturbation. The SF, 
in contrast, must be viewed as a genuine physical effect (even if a delicate one, as the 
SF too is gauge dependent — see below). There indeed exists an interpretation of the 
motion — we will discuss it later — wherein the particle moves freely on a geodesic of a 
certain smooth, perturbed spacetime, subject to no SF. However, in this description 
the smooth geometry is not the physical spacetime of the background -I- particle system 
(the metric of this geometry is not a retarded solution of the linearized Einstein 
equations), and so this effective spacetime cannot be said to represent "physical 
reality" any better than the SF itself. 

Thanks to work commencing around 1997 and continuing to these days, we now 
have a rather satisfactory theory of the gravitational SF in curved spacetime, which, 
in particular, addresses all of the above issues. Initial work was strongly inspired 
by the classical analyses of the analogous electromagnetic self-force problem in flat 
(Dirac, 1938 [6l]) and curved (DeWitt and Brehme, 1960 [58]) spacetimes. DeWitt 
and Brehme's work, especially, provided much of the mathematical framework — that 
of covariant bi-tensors in curved spacetime — needed for analyzing the gravitational 
problem too. In 1997 two independent groups published three independent derivations 
of the gravitational SF. Quinn and Wald jl41j used an axiomatic approach, in which 
the physical self-acceleration is deduced, essentially, by comparing the (divergent) self- 
field of the particle in question with that of a particle in a suitably constructed tangent 
flat space. Independently, Mino, Sasaki and Tanaka (MST) jll5j obtained the same 
result using a local energy-momentum conservation argument — a direct application 
of DeWitt and Brehme's method to the gravitational case. Both these methods pre- 
assume a notion of a point mass particle, and neither seeks to make a consistent sense 
of this notion. However, MST's paper also reported a second, independent derivation 
of the SF, using an approach which, for the first time, offered a fully GR-consistent 
treatment of the problem. 

This approach — a new implementation of the old idea of matched asymptotic 
expansions — relies on the assumption that there can be identified two separate 
lengthscales in the problem: One associated with the particle's mass and another, 
much larger, associated with the typical radius of curvature of the geometry in which 
the particle is moving. In the strong-field EMRI problem, the second lengthscale is 
provided by the mass of the central black hole, M ^ /x. MST's construction further 
assumes that the "particle" is actually a Schwarzschild black hole of mass /i. The 
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two separate scales in the setup define a "near zone", r <C Af, and a "far zone", 
r 3> ?n (wliere r is a suitable measure of distance from the small hole). In the near 
zone, the geometry is approximately that of the small Schwarzschild hole, with small 
tidal-type corrections from the background geometry. As we zoom away from the 
small object and enter the far zone, the effect of the small object's detailed structure 
becomes gradually less important, and at the far zone limit the geometry becomes 
that of the background spacetime, weakly perturbed by what is now a distant "point 
particle" — it is indeed the far-zone limit through which a notion of point mass can be 
defined in a consistent way. In situations where M ^ m one would have a "buffer 
zone" where m <^ r <^ M and both "near zone" and "far zone" descriptions of the 
geometry are valid. MST showed that matching the near zone and far zone metrics 
(expressed as asymptotic expansions in r/M and m/r, respectively) constrains the 
motion of the particle (from a far-zone point of view) and thus yields an expression 
for the SF. This expression agreed with those obtained by Quinn and Wald and by 
using DeWitt and Brehme's method; it was later coined the MiSaTaQuWa formula, 
an acronym based on the names of the five contributing authors. A self-contained and 
pedagogical review of the MiSaTaQuWa formula, including an elegant reproduction 
of previous derivations, can be found in Poisson's [131j . 

More recently, Gralla and Wald [77] (see also [78] for a more concise presentation) 
developed a new procedure for deriving the gravitational SF, which offers improved 
mathematical rigor as well as a generalization of the MiSaTaQuWa formula. Rather 
than relying on two separate asymptotic expansions of the metric as in MST's original 
method, Gralla and Wald introduce a single one-parameter family of metrics, which, 
through two different limiting procedures, can produce both near and far zone metrics 
in a natural way. This allows to define more robustly the criteria for existence of 
the two zones, and enables a more elegant buffer-zone matching. The analysis proves 
that in the far-zone limit the particle is described precisely by the usual delta- function 
distribution, and that at the very limit fj, ^ this particle moves on a geodesic of the 
background. Furthermore, the analysis relaxes all assumptions about the nature of 
the small object: It no longer need to be a Schwarzschild black hole, but can assume 
the form of any sufficiently small black hole or a blob of ordinary matter. This allows, 
in particular, for a spin-force term to appear in the resulting, generalized version of 
the MiSaTaQuWa formula. 

The main end product of the above theoretical developments is a firmly- 
established general formula for the SF in a class of background spacetimes including 
Kerr. It should be stressed that the SF formula stems in a deterministic way from 
nothing else than the Einstein equations with the usual conservation laws; it does 
not rely — and should not rely, as a matter of principle — on any form of ambiguous 
"regularization" or "subtraction of infinities" . This idea is implemented explicitly in 
the matched asymptotic expansions approach, and even more so in the Gralla- Wald 
analysis. 

2.3. MiSaTaQuWa equation 

We now state the MiSaTaQuWa formula [MTI dUl [m]— Eq. ^ below. (For 
simplicity we ignore spin-force terms and focus on the self-interaction part of the 
formula.) This will serve as a starting point for the rest of this review. 

Consider a timelike geodesic F in a background spacetime with metric gap. For 
concreteness, let us think of F as a test particle orbit outside a Kerr black hole, so gap 
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is the Kerr metric. Let r be the proper time along F, and let x" = -z"(t) describe F 
in some smooth coordinate system and u°' = dz^/dr be the four velocity of the test 
particle. Denote by /i™^ the physical, retarded metric perturbation from a particle of 
mass /X whose worldline is F. Assume is given in the Lorenz gauge: 

Vf'hTp = 0, (14) 

where /i™^, recall, is the trace-reversed version of /i^'^J [see Eq. ([7])]. Remember that 
throughout our discussion indices are raised and lowered using the background metric 
gafii and covariant derivatives are taken with respect to that metric. 

At any specetime point x, the retarded perturbation can be written as a sum of 
two pieces, 

Ki = h% + h>:i (15) 

the former being the "direct" contribution coming from the intersection of the past 
light-cone of x with F, and the latter being the "tail" contribution arising from the part 
of F inside this light cone (see Fig. [T|) . The occurrence of a tail term is a well-known 
feature of the wave equation in S-flD curved spacetime, and it can be interpreted 
physically as arising from the effect of waves being scattered off spacetime curvature 
("failure of the Huygens principle"). Both li^^p and h'^^ obviously diverge when 
evaluated on F; however, h}^^ is continuous and differentiable everywhere, including 
on the worldline. Notably, though, the tail field is not a smooth function on the 
worldline, and is not a vacuum solution of the linearized Einstein equations. 




Figure 1. An illustration of the setup described in the text. z{t) is a point on 
the timelike worldline F (thick solid line) and a; is a field point close to z, shown 
with a portion of its past light cone, e is the spatial geodesic distance from x to 
r and 5x°' = x°' — z". The metric perturbation at x consists of a direct and a 
tail contributions, illustrated by the thick dashed and dash-dot lines, respectively. 
[Graphics reproduced from Ref. I12j .] 



The MiSaTaQuWa formula states that the gravitational SF at a given point z 
along F results simply from the back reaction of the tail field: 

FS,u{z) = fiW'^^^fiz). (16) 

Here V"^''' is the usual "force operator", given in equation ([§]), which, recall, is 
dependent upon the four-velocity u" and the background metric ga/^ at point z||| 

t In the original MiSaTaQuWa formulation, the SF is not expressed directly in terms of the gradient 
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2.4- Detweiler- Whiting reformulation 

In a 2003 paper [S5| (see also a preliminary discussion in [52]) Detweiler and 
Whiting (DW) proposed an alternative formulation, which offers an interesting re- 
interpretation of the perturbed motion. DW replaced the direct/tail decomposition 
of the retarded perturbation [Eq. (fT5|) ] with a new decomposition, 

hT0 = hl(3 + h%, (17) 

where the i?-field h^^ (R for Regular), unlike the tail field, is a certain smooth, vacuum 
solution of the perturbation equations, which, nonetheless, gives rise to the same 
physical SF as the tail field: 

F,",if(z) = MV"^^/i^,(z). (18) 

The S'-field /i^^ (S for Singular), which precisely mimics the singular behavior of 
the retarded field near the particle, exerts no SF and does not affect the motion of 
the particle. The precise formal prescription for constructing the R and S fields is 
described nicely in Poisson's review |131| . 

DWs discovery that the MiSaTaQuWa SF can also be expressed as the back- 
reaction force from a smooth vacuum perturbation leads to an interesting re- 
interpretation of the gravitational SF effect: The particle effectively moves freely along 
a geodesic of a smooth perturbed spacetime with metric gap + h^p. In this alternative 
picture — more in the spirit of GR's equivalence principle — the notion of a SF becomes 
artificial (and obsolete) in much the the same way that the notion of an external 
gravitational force is artificial. 

It should be understood, however, that the i?-field does not represent the actual 
physical perturbation from the particle (the physical perturbation field is of course 
/i™J). The _R- field has peculiar causal properties which make it problematic as a 
candidate for what we may call a "physical field" : The value of the i?-field at an event 
X depends not only on events in the the causal past of x but also on events outside the 
light-cone of x |131j . Rather than an entity of physical substance, the i?- field should 
be viewed as an effective field that allows us to describe the dynamics in terms of 
geodesic motion. 

The two descriptions of the perturbed motion — self-accelerated motion in g vs. 
geodesic motion in g -f — are alternative (equivalent) interpretations of the same 
(genuine) physical effect. The two points of view are not contradictory but rather they 
are complementary in their perspective on the problem. Workers in the field often find 
it useful to invoke both descriptions alternately in order to get a fuller picture of the 
physics in question. Sago et al. [147J recently demonstrated the equivalence of the two 
approaches with an explicit calculation of a certain gauge-invariant physical SF effect 
in a particular example; we shall return to discuss this work in Sec. [H 

2. 5. Singular field 

Equations (fT6| and ([T8|) prescribe the correct regularization of the gravitational SF 
and form the fundamental basis for all modern SF calculations. In later sections we 

of /i^'^', but rather as an integral over the gradient of the relevant retarded Green's function along 
the portion of the worldline to the past of z [cf. Eq. (1.9.6) of Poisson 11311 ]. The commutation 
of the derivative and the worldline integral produces local terms at z, which, however, vanish upon 
contraction with V"^''. This leads to the equivalent formulation shown here in Eq. (I16I I. 
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will discuss practical methods for implementing these formulas. We will then need 
some more information on the properties of the direct (or singular) field, which we 
now give. 

The fields h'^'^ and share the same leading-order singular structure near the 
particle's worldline F, but they differ in their sub-dominant singular behavior. More 
precisely (referring again to Fig. [T]), consider a particular point z on F and a nearby 
off-F field point x. The form of both the direct and S fields is then given, in the 
Lorenz gauge, by [113 EH [121] 

TSAir/ ^ _ ^^iUa{x■,z)up{x■,z) ^ fJ- ^If" {x ; z) ^ g^^.^ 

7{^) + e{x; z) + "^"^ ' ^^^^ 

where Uf^ is the four-velocity vector parallelly-propagated from z to x, e is the spatial 
geodesic distance from a; to F (i.e., the proper length of the short normal geodesic 
section connecting x to F), w^'l^" are smooth functions of x (and z) which vanish 

at a; ^ z at least quadratically in the coordinate differences — z", and c^'^"^ are 
constants (dependent on z but not on x). The S and direct fields differ only in the 
explicit form of Wafj and (possibly) in the value of Cap, but neither of these will be 
important for us in this review. 

It is, however, important to emphasize that the singular form (jl9p is generally 
gauge dependent. The expression given here is specific to the Lorenz gauge, and is 
not guaranteed to keep its form if one makes other gauge choices for the perturbation. 
Indeed, it has been demonstrated [21] that some common gauge choices can "distort" 
the local rest-frame isotropy of the Lorenz-gauge singularity manifest (at leading order) 
in Eq. (fTg]) . 

2.6. Gauge dependence 

Earlier we emphasized some important differences between the concepts of a general 
(external) gravitational force and that of the SF. The two notions, however, share a 
basic common feature: They are both defined via a mapping of the physical trajectory 
from a "perturbed" spacetime to a "background" spacetime. In both cases, such 
a mapping procedure gives rise to a gauge ambiguity. A thorough analysis of the 
gauge dependence of the gravitational SF was presented in Ref. [H], and a gauge 
transformation law for the SF was derived. Consider again the infinitesimal gauge 
transformation of Eq. ([10|) . where now the gauge displacement vector ^" is assumed 
to scale like the particle's mass fx. The change this induces on the physical SF was 
found [21] to be given by 

= -M [(5"^ + n^u^)'(x + R'^.x.u^^ew] . (20) 

This has the same form as the gauge transformation law for the external gravitational 
force [recall Eq. (|13|) ]. which is not surprising given the similar geometrical origin of 
the gauge freedom in both cases. 

We note, however, that the derivation of the transformation rule ([20]) in Ref. [21] 
(which pre-dates the DW analysis) could not — and did not — follow the procedure we 
implemented above in deriving Eq. (|13p . simply replacing the external perturbation 
hap with the physical perturbation from the particle, /i™^. That is because /i™^, unlike 
the smooth external field hafi, is singular at the particle, and so expressions such as 
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(jlip or (fT^ (with /iq^ ^ ^^a/s) would make no sense when evaluated along the particle's 
worldline. DW's later i?-field interpretation suggests an alternative derivation of Eq. 
(|20p: Since the i?- field can be viewed, effectively, as a smooth external perturbation 
(even if one with peculiar causal properties), and since the SF is the force exerted 
by this perturbation, the derivation leading to Eq. (fT3l) can be repeated in full with 
hai3 /ij^^, and Eq. ([20|) follows immediately. In effect, the i?- field interpretation 
allows us here to treat the SF on an equal footing with the gravitational force from 
an external perturbation. 

The gauge dependence of the SF by no means implies that there is something 
"unphysical" about it — the SF is as physical as the metric perturbation itself, which 
is also gauge dependent. The gauge dependence does mean, however, that one needs 
to exercise some care in decoding the physical content of the SF. One cannot expect 
to be able to describe the physical effect of the SF based on the value of the SF alone 
(to put this to extreme: one can always make a gauge choice that nullifies the SF 
anywhere along the orbit!). Instead, a meaningful description of the physical effect 
must involve both the SF and the gauge information associated with it (in the form 
of the metric perturbation, for example). 

Another crucial point to have in mind is that the MiSaTaQuWa formula (fT6|l is 
guaranteed to hold true only if h'^^p satisfies the Lorenz-gauge condition p4|l . Strictly 
speaking, the SF in a given non-Lorenz gauge only makes sense if, for relating that 
gauge to the Lorenz gauge, the expression on the right-hand side of Eq. ([20]) is well 
defined. This is not at all an obvious condition; some simple counter-examples are 
analyzed in [21] . In gauges for which the right-hand side of ((20|) does not have a well 
defined (i.e., finite and direction independent) particle limit, one might still devise a 
useful notion of the SF by averaging over angular directions, or by taking a directional 
limit in a consistent fashion. Such possibilities are discussed in Refs. [2T1 177] . 



2. 7. Equations of motion 

Given the SF, the particle's equation of motion becomes 

^iu''Vpu'^ = F,%, (21) 

where on the left-hand side we have the usual four-acceleration (times /i) along the 
background trajectory. This equation, along with the MiSaTaQuWa equation (fTB]) . 
describe the dynamics of the particle given the metric perturbation (and assuming 
one has a way of extracting the tail piece out of the full perturbation). To close 
the system of equations we need to know how the metric perturbation is determined 
from the particle's trajectory. This, of course, is provided by the linearized Einstein 
equation, which takes the Lorenz-gauge form 

/oo 
(_5)-i/2 j4^^,, _ dr, (22) 

-oo 

where g is the determinant oi gafs and z^{t), recall, describes the particle's worldline. 
The source on the right-hand side is the usual distributional representation of the 
point particle's energy-momentum. The field equation (|22p is to be supplemented by 
the gauge condition p4p and by suitable boundary conditions. 

The set of equations p6|) , ((2T|) , ((22|) and p4|) (together with a method to obtain 
^taii q£ h'^^^) should in principle determine the dynamics of the orbit at linear order 
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in perturbation theory. However, as pointed out by Garlla and Wald recently |77j . 
the field equation (|^^ is only consistent with the Lorenz gauge condition (fHl) if the 
particle is moving strictly along a geodesic — which would then be inconsistent with 
the equation of motion (|2ip . To resolve this inconsistency while allowing for orbital 
evolution, Gralla and Wald suggested a "Lorenz-gauge relaxation" approach, wherein 
one relaxes the gauge condition p4|l and considers solutions of the set (|16I21I22|) . 
One then expects that, in situations where the orbit is very nearly geodesic (as is 
usually the case with LISA-relevant astrophysical inspirals) , such solutions would give 
a faithful, albeit approximate description of the actual orbit. This approach is yet to 
implemented and tested in actual calculations of the orbital evolution. 

The proposal to use self-consistent solutions of the set (|16[ I21I22P for modeling 
the slow orbital evolution at linear order in perturbation theory was put forward by 
Gralla and Wald in [77]. A different mathematical framework, based on techniques 
from multi-scale perturbation theory, was developed by Hinderer and Flanagan in 
[57] . Pound and Poisson |136l I137j performed first actual calculations of the orbital 
evolution under the full SF effect, using multi-scale analysis, within a weak-field, PN 
framework. Thus far there are no calculations of the full orbital evolution in strong- 
field scenarios. 

In the rest of this review we will not consider any further the question of orbital 
evolution, but rather focus on the calculation of the SF [via Eqs. (HHJ) or p^ ] along a 
pre-determined, non-evolving orbit. We consider this a preliminary step, whose output 
(e.g., the value of the SF at sufficiently many points across the relevant phase space) 
can later be incorporated into whichever evolution scheme one chooses to apply. 



3. Overview of implementation framew^orks and calculations to date 

Starting in the late 1990s, work began to translate MiSaTaQuWa's SF formalism 
into practical working schemes and to implement it in actual calculations. While 
the "holy grail" of this program has from the outset been — and still remains — the 
calculation of the gravitational SF for generic orbits around Kerr black holes, much 
of the initial effort has focused on the simpler toy problem of the scalar-field SF, 
and on simple classes of orbits (radial, circular) in Schwarzschild spacetime. The last 
few years, however, saw first calculations of the gravitational (and electromagnetic) 
SFs for generic orbits in Schwarzschild — and initial work on Kerr is now underway. 
In this section we give a broad overview of the methods that have been suggested 
for implementing MiSaTaQuWa's formalism, and we survey the actual calculations 
performed so far. 

Our starting point are the formal expressions ([16]) and (jlSp for the gravitational 
SF. To facilitate the following discussion, we first use Eqs. ([T5| and p7ll to recast 
these expressions in the more practical form 

FS,u{z) = M lim [^f-^hl'^ix] - K'^'hffi'^)] ' (23) 

referring once again to Fig.[l]for notation. This describes a "regularization" procedure, 
which one can perform using either ft,^"' or h^^ — both producing the same final value 
for the SF. The limit procedure is necessary here because the individual fields h^^^ 
and Ti^"^^ , unlike their difference Tij^^^^^, are each singular at a; — > z, and so are their 
derivatives. Since the operator V"^''' is only defined along the particle's worldline 
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[as it involves the four- velocity — recall Eq. ([5])], in Eq. we needed to introduce 
an extension of this operator off the worldline — denoted V"^'''. For all x in the 
neighborhood of a given worldline point z, the operator is given by Eq. 

where g"^ and u" take the same values they have at z, and Vs is the standard 
covariant derivative at point x. Here, and throughout this review, we use this "fixed 
contravariant components" extension of V"''''' exclusively, although other natural 
extensions are possible [23]. Of course, the final value of the SF in Eq. (|23p. after 
the limit x — > z is taken, is not sensitive to the choice of extension. Notice that 
the definition of V"''''' is, of course, coordinate dependent, and it becomes well-posed 
only in reference to a specific coordinate system. Also note that the operator W"^'^ is 
defined with respect to a given point z. 

The most basic technical challenge one faces in preparing to implement the 
MiSaTaQuWa formula is the so-called subtraction problem: How does one go about 
extracting the tail (or R) piece from the full retarded perturbation in practice? 
Equation (|23p suggests applying the subtraction /i™^ — ti^^^^, using approximate 
analytic expressions for the direct/S fields, such as the one in Eq. p9|) . However, 
this involves the removal of one divergent quantity from another, which is not easily 
tractable in actual numerical calculations. Several strategies have been proposed to 
address this problem, and in the main part of this review we shall describe a few of 
them in some detail. Here we proceed with a brief overview of the main avenues of 
approach to this problem. 

Quasi-local/matched expansions calculations. [H [71 [6l [5l I126[ 1127] This approach 
tackles the calculation of the tail contribution directly, by analytically evaluating the 
Hadamard expansion of the relevant Green's function. Such calculations capture 
the "near" part of the tail, which, one might hope, represents the dominant 
contribution in problems of interest. Quasi- local calculations can be supplemented 
by a numerical computation of the "far" part of the tail, a strategy referred to as 
"matched expansions" (not to be confused with matched asymptotic expansions) . The 
applicability of this idea was demonstrated very recently [48] with a full calculation 
in Nariai spacetime (a simple toy spacetime featuring many of the characteristics of 
Schwarzschild) . 

Weak-field analysis. [SSI 11291 11351 11361 I137j The tail formula can be evaluated 
analytically for certain weak-field configurations, within a Newtonian or a post- 
Newtonian (PN) framework. Such work has provided important insight into the nature 
and properties of the SF. PN techniques have also been implemented in combination 
with the mode-sum method discussed below [III [Ml [551155]. 

Radiation-gauge regularization. [96| This approach proposes a reformulation of the 
MiSaTaQuWa regularization, in which one reconstructs the R-part of the metric 
perturbation in a radiation gauge (rather than in the Lorenz gauge) from a suitably 
regularized Newman-Penrose curvature scalar (ipo or ipi)- The main advantage of this 
method is that it reduces the numerical component of the calculation to a solution of 
a single scalar-like (Teukolsky's) equation. However, some of the technical complexity 
is relegated to the metric reconstruction step. This technique was implemented so far 
only for a particle held static in Schwarzschild, but more interesting cases are currently 
being studied [95l[67]. See the discussion in Sec. 15.11 for more details. 
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Mode-sum method. [H \IM [H El [HISlllllMllSTllHSllHilMl An approach 
whereby one evaluates the tail contribution mode by mode in a multipole expansion. 
The subtraction "ret— dir" (or "ret— S") in Eq. ([^5)) is performed mode by mode, 
avoiding the need to deal with divergent quantities. The method exploits the 
separability of the field equations in Kerr into multipole harmonics. The mode-sum 
method has provided the framework for the bulk of work on SF calculations over the 
last decade. We will discuss it in detail in Sec. [H 

Puncture methods, p!? ! [28 l fT58l [95l I105j A set of recently-proposed methods custom- 
built for time-domain numerical implementation in 2+1 or 3+1 dimensions. Common 
to these methods is the idea to utilize as a variable for the numerical time-evolution 
a "punctured" field, constructed from the full (retarded) field by removing a suitable 
singular piece, given analytically. The piece removed approximates the correct S-field 
sufficiently well that the resulting "residual" field is guaranteed to yield the correct 
MiSaTaQuWa SF. In the 2+lD version of this approach the regularization is done 
mode by mode in the azimuthal (m-mode) expansion of the full field. This procedure 
offers significant simplification; we shall review it in detail in Sec. [71 

3.1. SF calculations to date 

As we have mentioned already, the program to calculate the SF for black hole orbits has 
been progressing gradually, through the study of a set of simplified model problems. 
Some of the necessary computational techniques were first tested within the simpler 
framework of a scalar-field toy model before being applied to the electromagnetic (EM) 
and gravitational problems. Authors have considered special classes of orbits (static, 
radial, circular) before attempting more generic cases, and much of the work so far 
has focused on Schwarzschild orbits. The state of the art in the field are numerical 
codes to compute the scalar, EM and gravitational SFs along any given (geodesic) 
orbit outside a Schwarzschild black hole. It is reasonable to expect that attention will 
now be increasingly drawn to the Kerr problem. 

The information in Tables [TH3] is meant to provide a quick reference to work done 
so far. It covers actual evaluations of the local SF that are based on the MiSaTaQuWa 
formulation (or the analogous scalar-field and EM formulations of Refs. |.140j and 
[58l [88l I141j . respectively), either directly or through one of the aforementioned 
implementation methods. We have included weak-field and FN implementations, but 
have not included work based on the radiative field approach. Some of the techniques 
referred to under "strategy" will be discussed in the following sections. 

4. Mode-sum method 

Let us write the subtraction formula (|23p using the more compact notation 



0(z)= hm [FZdx)~F^{x)], 



(24) 




(25) 



For concreteness and simplicity we adopt here the S-field subtraction, noting that the 
entire discussion in this section would not be altered upon replacing S^direct and 
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Scalar-field self force 


Case 


Author (s) 


Strategy 


Newtonian potential 


Pfennine & Poisson [129] 


direct, analytic 


spherical mass shell: 


Burko et al. [15] 


mode sum, analytic 


static particle 






isotropic cosmology: 


Burko et al. [44] 


direct, analytic 


static particle 






isotropic cosmology: 


Haas & Poisson [82j 


direct, analytic 


slow motion 






Nariai spacetime: 


Casals et al. [35] 


matched expansions 


static particle 






Schwarzschild: 


Burko [39 


mode sum, analytic 


static particle 


Wiseman [163 


direct, analytic 


Schwarzschild: 


Barack & Burko [13] 


mode sum, numerical 


radial geodesies 




(I+ID evolution) 


Schwarzschild: 


Nakano et al. ]T5D] 


post-Newtonian, analytic 


circular geodesies 


Hikida et al. [86] 






Burko [40] 

! 1 


mode sum, numerical 




Detweiler et al. [57, 60 


ffrequencv domain! 




Vega & Detweiler [158] 


puncture, numerical 
(1+lD evolution) 




vega et al. [159] 


puncture, numerical 
(3+ ID evolution) 


Schwarzschild: 


Haas [HO] 


mode sum, numerical 


eccentric geodesies 




(1+lD evolution) 


Kerr-Newman: 


Burko & Liu [42] 


mode sum, analytic 


static particle 






Kerr: circular- 


Warburton & Barack [162( 


mode sum, numerical 


equatorial geodesies 




(frequency domain) 



Table 1. Calculations of the scalar-field SF as a toy model for the gravitational 
problem. In this table, as well as in Tables [2] and \3\ "direct" implies explicit 
evaluation of the tail contribution to the SF. "Mode sum" and "puncture" refer 
to the two implementation schemes described, respectively, in Sec. |4] and [7] of 
this review. The method of "matched expansions" (different from "matched 
asymptotic expansions") is described briefly in Secs.|3] 



R^taili The fields F^^^ and Fg inherit the extension ambiguity of V"'''*', but here 
we shall always use the (coordinate dependent) "fixed" extension described above. 
Both F^^^ and Fg, of course, diverge at the particle, but their difference is a smooth 
(analytic) function of x even at the particle. 

In the mode-sum method we formally decompose each vectorial component 
of both F^^j. and F^ into spherical harmonics. These harmonics are defined in 
the Kerr/ Schwarzschild background based on the Boyer-Lindquist/Schwarzschild 
coordinates {t,r,9,ip) in the standard way, i.e., through a projection onto an 

I The only exception is that statements referring to the smoothness of the R field would need to be 
formulated more carefully to reflect the irregularity in the higher derivatives of the tail field. This 
irregularity, however, would have little practical impact on the discussion in this section. 
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Electromagnetic self force 


Case 


Author (s) 


Strategy 


Newtonian potential 


Pfenning & Poisson 129J 


direct, analytic 


Isotropic cosmology: 


Haas & Poisson [82] 


direct, analytic 


slow motion 






Schwarzschild: 


Burko [391 


mode sum, analytic 


static particle 








Keidl et al. 96^ 


radiation-gauge 
regularization, analytic 


Schwarzschild: 


Haas [ST] 


mode sum, numerical 


eccentric geodesies 




(1+lD evolution) 



Table 2. Recent calculations of the Electromagnetic self-force. 



Gravitational self force 


Case 


Author (s) 


Strategy 


Newtonian potential 


Pfenning & Poisson [129] 


direct, analytic 


Schwarzschild: 
radial geodesies 


Barack & Lousto [18j 


mode sum, 1-1- ID evolution 
in Regge- Wheeler gauge 


Schwarzschild: 
static particle 


Keidl et al. [96] 


radiation-gauge 
regularization, analytic 


Schwarzschild: 
circular geodesies 


Barack & Sago [25] 


mode sum, 1-f ID evolution 
in Lorenz gauge 




Dctweiler [54 


mode sum, frequency-domain 
in Regge- Wheeler gauge 


Schwarzschild: 
eccentric geodesies 


Barack & Sago [26] [143] 


mode sum, 1-1- ID evolution 
in Lorenz gauge 



Table 3. Calculations of the Gravitational self-force. 



orthogonal basis of angular functions defined on surfaces of constant t and r. Let 
us denote by F^^^ {x) and F^' [x) the Z-mode contributions to F^"^ and , respectively 
(summed over m). A key observation is that each of these Z-mode fields is finite even 
at the particle. This suggests a natural regularization procedure, which, essentially, 
amounts to performing the subtraction F^"^ — F^ mode by mode. The idea is best 
developed through an elementary example, as follows. 

4.1. An elementary example 

Consider a pointlike particle of mass fi at rest in flat space. The location of the 
particle is x = Xp in a given Cartesian system. In this simple static configuration the 
perturbed Einstein equations (|22p read 



— 167r/i S^(x — Xp), 



(26) 
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where is the 3D Laplacian, and with all other components of h}^^ vanishingl§| 
The static perturbation h^^p automatically satisfies the Lorenz-gauge condition (fH)) . 
Of course, in this simple case we can immediately write down the exact physical 
(Coulomb-like) static solution, h^p = A^/\x — Xp\, and we also trivially have F^^f = 0. 
However, for the sake of our discussion, let us proceed by considering the multipole 
expansion of the perturbation. 

To this end, introduce polar coordinates (r, 9, ip), such that our particle is located 
at Xp = (ro 7^ 0, Vo)) f^nd expand the physical solution hn in spherical harmonics 
on the spheres r =:const, in the form 



1^0 in~~l 



h\f ^y-h[,{r,e), where h\,{T,e) = y~h[r{r)YUO,v)- (27) 



This expansion separates the field equation ([26| into radial an angular parts, the 
former reading (for each /, m) 

h\Z, + -h\Z - h\r = YC^{e,, ^o) Sir - ro), (28) 

r r 

where an asterisk denotes complex conjugation and a comma denotes partial 
differentiation. The unique physical I, m-mode solution, continuous everywhere and 
regular at both r = and r ^ oo, reads 



givmg 



^"W = ^^>^4(^o,.o)x|;-;. ' (29) 



/.j,(r,^)^i^p.(cos7)x| (30 



(r/ro)', r < ro. 



where Pi is the Legendre polynomial and 7 is the angle subtended by the two radius- 
vectors to X and Xp (see Fig. [5]). 

We now construct the force field F^^^ as it is defined in Eq. (1^ . We find = 0, 
and the spatial components are F^^^ = fjS/'^^^'htt = {iJ,/4:)h^f. Focus now on the r 
component. The I mode of F^^^ is given simply as F^^^ = {ii/A)h[^ ^. Using Eq. ([SO]) 
and evaluating F^'^j at the particle (taking 7^0 followed by r — > r^), we obtain 

F;,',±(fp) = Ti4-#i' ^here L^l + \. (31) 

Here the subscripts ± indicate the two (different) values obtained by taking the particle 
limit from "outside" (r ^ r"*") and "inside" (r r~). 

Let us note the following features manifest in the above simple analysis: 

• The individual I modes of the metric perturbation, ft.^^, are each continuous at 
the particle's location, although their derivatives are discontinuous there. 

• The individual I modes -Fj-oti have finite one-sided values at the particle. 

§ Here the label 'ret' is less appropriate, but we shall retain it to adhere to our general notation. As 
in the general case, h^^p represents the unique physical perturbation — here the unique regular static 
solution of Eq. II26I I. 



L.Barack 



20 



■"■fi-*.Xp=(;-,e„,<p„) 




/ 




X = (r,9,(p) 



Figure 2. An illustration of the simple setup described in the text: A particle of 
mass fj, in flat space is at rest at Xp = {ro,Oo,ipo). The gravitational field of the 
particle is decomposed into spherical harmonics, each contributing a finite amount 
to the full radial force acting on the particle: either F^_i_ or F^_, depending on 
whether the force is calculated from r ^ r+ or r — > r~. x = {r,6,ip) is an 
arbitrary field point used in the construction described in the text. 



• At large I, each of the one-sided vahies of at the particle is dominated by 
a term cx I. (The mode sum obviously diverges at the particle, reflecting the 
divergence of the full force i^^ot there.) 

It turns out (as we shall see later) that all above features are quite generic, and 
they carry over intact to the much more general problem of a particle moving in Kerr 
spacetime. Specifically, we find that, at any point along the particle's trajectory, the 
(one-sided values of the) modes always admit the laige-l form 



In our elementary problem the power series in 1/L truncates at the term, but 
in general the series can be infinite. The ^-independent coefficients Aa, and Cq., 
whose values depend on the background geometry as well on the particle's location 
and four- velocity, are characteristic of the local structure of the particle singularity at 
large I. These coefficients, called Regularization Parameters, play a crucial role in the 
mode-sum regularization procedure, as we describe next. 

4.2. The mode- sum formula 

Consider a mass particle moving on a geodesic trajectory in Kerr, and suppose we 
are interested in the value of the SF at a point z along the trajectory, with Boyer- 
Lindquist coordinates (to, fo, ^0, <^o)- Starting with Eq. (j^ . let us formally expand 
F^^^(x) and F^{x) in spherical harmonics on the surfaces t, r=const. Here we ignore 
the vectorial nature of i^^gj and F^ and, for mathematical simplicity, treat each of their 
Boyer-Lindquist components as a scalar function (see [53] for a more sophisticated, 
covariant treatment). Denoting the respective Z-mode contributions (summed over m) 
by F^J;^{x) and F^\x), we write 



F' 



rctib ~ 



(32) 



00 




(33) 
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We remind that the form of F^i and will depend on the specific off-worldhne 
extension chosen for V"''''; however, this ambiguity disappears upon taking the hmit 
X ^ z, and the final SF is, of course, insensitive to the choice of extension. 

Since F^^^{x) — F^{x) is a smooth function for all x, the mode sum in Eq. ([55]) is 
guaranteed to converge exponentially for all x (this is a general mathematical property 
of the multipole expansion). In particular, the sum converges uniformly at a; = z, and 
we are allowed to change the order of limit and summation. We expect, however, that 
the particle limit of the individual terms F^^. and Fg*' is only defined in a one-sided 
sense, as in the elementary example studied above. Hence, we write 

oc 

^self (^) = E [^rZL (^) - FS± (^)] , (34) 
1=0 

where ± indicates the values obtained by first taking the limits t ^ to, 9 ^ Oq and 
tp ^ Lpo, and then taking r ^ . Of course, the difference F^[ _(_(z) — Fg*j_(z) does 
not depend on the direction from which the radial limit is taken: As the difference 
F^{{x) — Fg\x) is the / mode of a smooth function, it is itself smooth for all x. 

Furthermore, since the mode sum in Eq. (j34p converges faster than any power of 
1/L, we expect Fj."' and i^g*' to share the same large-Z power expansion ([5^ . with the 
same expansion coefficients. This motivates us to re-express Eq. ([34]) in the form 

oo 

O(^) = E [^"et±(^) T LA- C^/L] 

oc 

- E [^s"± (^) T LA^ - - CyL] . (35) 

1=0 

Here, each of the terms in square brackets falls off at least as ~ as I oo, and 
so each of the two sums converges at least as ^ 1/L We hence arrive at the following 
mode-sum reformulation of the MiSaTaQuWa equation: 

oo 

^^scif(^) = E Kcl±(^) T LA" C"] - D", (36) 

1=0 

with 

oo 

£»" = E [^s±(-2) T LA" -B"- C"] . (37) 

(=0 

The mode-sum formula (jH]), first proposed in Refs. [5D] (scalar-field case) and 
[TT] (gravitational case) provides a practical way of calculating the SF, once the 
regularization parameters A" , B" , C" and D" are known. It is practical because (i) it 
involves no subtraction of divergent quantities, and (ii) it builds naturally on the fact 
that in black hole perturbation theory one usually calculates the metric perturbation 
mode by mode in a multipole decomposition. 

The values of the regularization parameters are obtained analytically via a local 
analysis of the singular (or direct) field near the particle. Early calculations of these 
parameters [20l [ini E] , which were restricted to specific orbits, were based on a local 
analysis of the Green's function near coincidence {x — > z) at large Later, after the 
local form of the direct field had been derived explicitly to sufficient accuracy |116j . 
workers have been able to derive general expressions for the regularization parameters, 
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valid for generic (geodesic) orbits — first in Schwarzschild [571 [HI i then in Kerr |24j . 
Later improvements (in the scalar case) are due to Ref. [ST] , where higher-order terms 
in the 1 /L expansion were derived analytically in order to accelerate the convergence 
of the mode sum; and Ref. [83] , where the regularization parameters were redefined 
as scalar quantities using a covariant-form projection of the singular field onto a null 
tetrad based on the worldline. 

With the regularization parameters given in analytic form, the calculation of the 
SF within the mode-sum scheme follows this procedure: 

(i) For a given geodesic orbit, calculate sufficiently many multipole modes of the 
physical, retarded metric perturbation in the Lorenz gauge. (Here "multipole 
modes" may refer to tensor harmonics, Fourier-tensor harmonics, or spheroidal 
harmonics, depending on the case considered and on the perturbation framework 
used for tackling the field equations.) How this might be (and is being) done in 
practice — usually using numerical methods — will be discussed in Sec. [5] 

(ii) Construct the Z-modes F^[^ at the particle. An example of how this is done in 
practice will be provided in Sec. O 

(iii) Use the modes F^^j_{z) as input for the mode-sum formula p6|l . 

The mode-sum procedure has the very useful ability to automatically test itself: 
If either the values of the regularization parameters A", B" and C", or the values of 
the Z-modes -Fj-oti (usually computed numerically) are wrongly calculated, the l-mode 
sum in Eq. (|36|) will very likely fail to converge. Conversely, if one finds that the 
calculated i-term in the sum has the expected ~ I/i^ fall-off at large L, that provides 
an excellent validation test for the entire calculation. So far, the analytic values of 
the parameters A"', B"' and C" have been successfully "tested" (in the above sense) 
against numerical calculations for generic orbits in Schwarzschild, and also (in the 
scalar case) for circular equatorial orbits in Kerr. 

In what follows we state the values of the regularization parameters in the 
gravitational case, for generic geodesies in Kerr, as derived in Ref. [21]. In the 
Appendix we provide a full derivation of these parameters (not included in Ref. f24'). 

In passing, we remind that the mode-sum formula (in the gravitational case) is 
formulated in the Lorenz gauge, just like the MiSaTaQuWa formula on which it relies. 
The values of the regularization parameters; the form of the large-Z expansion in Eq. 
P2p : or even the very definiteness of the SF — all of these are gauge dependent. It 
has been shown, however, that the regularization parameters remain invariant under 
gauge transformations (from the Lorenz gauge) that are sufficiently regular (21| . 

4-3. Regularization parameters for generic orbits in Kerr 

We state here the values of the regularization parameters for the gravitational SF at an 
arbitrary point z along an arbitrary geodesic in Kerr spacetime. We assume the point 
z has Boyer-Lindquist coordinates (tQ,rQ,9o,'fio)- The regularization parameters C" 
and D°' are always zero: 

C" = £>" = 0. (38) 
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The Boyer-Lindquist components of the regularization parameter A°' at z are given 
by 

9 / -In \ 1/2 

A* = - = = 0, (39) 

where 

V = l + ul /gee + /g^^ . (40) 

Here is the particle's four-velocity and gai3 is the background Kerr metric — both 
at point z. 

The value of the parameter i?" is more complicated. It can be expressed in a 
compact form as 

B« = M'(27r)-ip«,,,7«''^^ (41) 

where hereafter reman indices (a, 6, c, ...) run over the two Boyer-Lindquist angular 
coordinates 0, (p. The cocfBcients P"abcd given by 

Plbcd = I [P"d{^Pabc + 2PabPc) - P'^H^Pxab + Pabx)Pcd] 

+ {SPlPbe - P%Pab)C'cd, (42) 

where 

Pa = U^uPgxp,a, Pap = gap + UaUp, Pafij = T^pPx^, (43) 

and 

= i sin ^0 cos ^0, ^^^o^-l cot ^o, (44) 

with all other coefficients C^^ vanishing, and with F^^ being the background 
connection coefficients at z. The quantities i"''^'' are 

r'^-'' = (sin^o)-^ r G(x)-'/'(sinx)^(cosx)'-^rfx, (45) 
Jo 

where 

G{x) = Pee cos^ X + '^Pe<p sin x cos x/ sin + P^^ sin^ x/ sin^ 6^0, (46) 

and = N{abcd) is the number of times the index tp occurs in the combination 
(a, b, c, d), namely 

TV = 5" +5^ + 5^ + 5^. (47) 

We may write /"^^'* explicitly in terms of standard complete Elliptic integrals. 
Introducing the short-hand notation 



a = sin^ 00 Pee/P^f - 1> /? = 2 sin 6»o Pe^/P^^,, (48) 

we have 

Qlf^k{w)-rlf^E{w) , (49) 



jabcd ^ (sin^o) ^ 



where 

Q = a + 2-(a2+/32)V2, (50) 
d = 3P'/?{smeo)-'{a^ + P^f{4a + 4 - l3^f/\Q/2y/^, (51) 
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K{w) = Jg^'^i^ — wsin^x) ^/'^dx and E{w) = — wsiv? x^/'^dx are complete 

Elliptic integrals of the 1st and 2nd kinds, respectively, and the argument is 

. . ^("^+^^)^^^ (52) 

The ten coefficients /jf ^ in Eq. ^ are given by 

/]^^ = 4 [I2a^ + a2(8 - 3/3^) - Aa[3^ + - 8)] , 

4")= -i6[8a3 + a2(4-7/32) + a/32(/32-4)-/32(/32 + 4)]^ (53) 

= 8/3 [9^2 - 2a(/32 - 4) + (5^] , 

^ - 4/3[12a3 - _ 52) + a(32 - 12/32) ^ p2^^p2 ^ 4^]^ (-54-) 

= _ 4 [8^3 - a2(/32 - 8) - 8aP^ + p^(3fi^ - 8)] , 

4^' = 8[4a4 + a^{/3^ + 12) - 2a^{P^ - 4) 

+ 3a/3' iP^ - 4) + 2/32 (3/32 - 4)] , (55) 

iP = 8/3 [a^ - 7a^ + a{3(3^ - 8) + p"^] , 

4) = - 4/3[8a^ - 4a3 + 0^(15/32 - 44) + Aa^bf]^ - 8) 

+ /32(3/32+4)], (56) 



4' = - 4[4a4 - 4a3 + 0^(7/32 - 8) + 12a/32 - /32(/32 - 8)], 
1^ = 16[4a^ + 4a4 + 0^(7/32 - 4) + a'^{lip^ - 4) 



(2a + l)/3"(/3"+4)]. (57) 



The sharp-eyed reader may observe that the expression for P^g^i,^^ in Eq. 
differs somewhat from that given in Ref. [24j . The reason for this discrepancy is 
that in Ref. [2^ we have made a different choice of extension for V"'''' [one in which 
the metric functions 

gap in Eq. ([9]) take their actual value at x, and is parallely- 
propagated from the worldline to x along a normal geodesic] . This affects the value of 
the parameter _B" only. We have opted here to give the parameter values corresponding 
to the "fixed" extension because that extension is more easily implemented in actual 
numerical calculations. To obtain the value of B" in the extension of Ref. [H], all one 
needs to do is replace the expression in Eq. ((42|) with 



PIbcAPP] - I [iP^dPabc - i2Pl, + Pab^Pcd] + {iPIPbe " PZPab)C^cd, (58) 

where 'PP' denotes the extension of Ref. [H] (which involves a Parallel Propagation 
of the four- velocity) . Interestingly, the regularization parameters for the scalar-field 
SF are precisely the same as the gravitational-field ones, with the parameter 
calculated using P1,bcd{PP] (^^"^ ^i^h the obvious replacement of the mass ^ with the 
scalar charge) [24]. 
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5. Numerical implementation strategies 

There are two broad (somewhat related) issues that need to be addressed in preparing 
to implement the mode-sum formula in practice (these issues arise, in some form, 
whether one uses the mode sum scheme or any of the other implementation methods 
based on the MiSaTaQuWa formalism). The first practical issue has often been 
referred to as the gauge problem: The mode-sum method (as the MiSaTaQuWa 
formalism underpinning it) is formulated in terms of the metric perturbation in 
the Lorenz gauge, while standard methods in black hole perturbation theory are 
formulated in other gauges. The second practical issue concerns the numerical 
treatment of the point-particle singularity; the problem takes a different form in the 
frequency domain and in the time domain, and we shall discuss these two frameworks 
separately below. We start, however, with a discussion of the gauge problem and the 
methods developed to address it. 

5.1. Overcoming the gauge problem 

The calculation of the gravitational SF requires direct information about the local 
metric perturbation near (and at) the particle. More specifically, one needs to be 
able to construct the metric perturbation, along with its derivatives, in the Lorenz 
gauge, at the particle's location. The mode-sum approach allows us to do so without 
encountering infinities by considering individual multipoles of the perturbation, but 
the problem remains how to obtain these multipoles in the desired Lorenz gauge. 
Unfortunately, standard formulations of black hole perturbations employ other gauges, 
favored for their algebraic simplicity. Such gauges are simple because they refiect 
faithfully the global symmetries of the underlying black hole spacetimes — unlike the 
Lorenz gauge, which is suitable for describing the locally-isotropic particle singularity, 
but complies less well with the global symmetry of the background. 

A common gauge choice for perturbation studies in Schwarzschild is the one 
introduced many years ago by Regge and Wheeler [142] . and further developed by 
Zerilli [164j and Moncrief [118j . In this gauge, certain projections of the metric 
perturbation onto a tensor-harmonic basis are taken to vanish, which results in 
a significant simplification of the perturbation equations. Another such useful 
"algebraic" gauge is the radiation gauge introduced by Chrzanowski |49j , in which 
one sets to zero the projection of the perturbation along a principle null direction 
of the background black hole geometry. Perturbations of the Kerr geometry have 
been studied almost exclusively using the powerful formalism by Teukolsky |155j . in 
which the perturbation is formulated in terms of the Newman-Penrose gauge-invariant 
scalars, rather than the metric. A reconstruction procedure for the metric perturbation 
out of the Teukolsky variables (in vacuum) was prescribed by Chrzanowski [49j (with 
later supplements by Wald |160j and Ori |123j ). but only in the radiation gauge. It is 
not known how to similarly reconstruct the metric in the Lorenz gauge. 

Several strategies have been proposed for dealing with this gauge-related difficulty. 
Some involve a deviation from the original MiSaTaQuWa notion of SF, while others 
seek to tackle the calculation of the Lorenz-gauge perturbation directly. Here is a 
survey of the main strategies. 

Self force in a "hybrid" gauge. Equation (I20p in Sec. [2] describes the gauge 
transformation of the SF. Let us refer here to a certain gauge as "regular" if the 
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transformation from the Lorenz gauge yields a well-defined SF in that gauge [this 
requires that the expression in square brackets in Eq. (|^D|) admits a definite and finite 
particle limit] . It has been shown [5T] that the mode-sum formula maintains its form 
(|36p. with the same regularization parameters, for any such regular gauges. Namely, 
for the SF in a specific regular gauge 'reg' we have the mode-sum formula 

oo 

^Tcif [reg] = [^rct± [reg] T LA^ - B^] , (59) 

where the force modes [reg] are these constructed from the retarded metric 

perturbation in the 'reg' gauge (we have already set C" = = here). 

Now recall the radiative inspiral problem which motivates us here: Although 
the momentary SF is gauge-dependent, the long-term radiative evolution of the orbit 
(as expressed, for example, through the drift of the constants of motion) has gauge- 
invariant characteristics that should be accessible from the SF in whatever regular 
gauge. And so, insofar as the physical inspiral problem is concerned, one might have 
hoped to circumvent the gauge problem by simply evaluating the SF in any gauge 
which is both regular and practical, using Eq. ([59)1 . Unfortunately, while the Lorenz 
gauge itself is "regular but not practical", both the Regge- Wheeler gauge and the 
radiation gauge are generally "practical but not regular", as demonstrated in Ref. 

However, a practical solution now suggests itself: Devise a gauge which is regular 
in the above sense, and yet practical in that it relates to one of the "practical" 
gauges — say, the radiation gauge — through a simple, explicit gauge transformation 
(unlike the Lorenz gauge itself). Heuristically, one may picture such a "hybrid" gauge 
(also referred to as "intermediate" gauge [24j) as one in which the metric perturbation 
retains its isotropic Lorenz-like form near the particle, while away from the particle it 
deforms so as to resemble the radiation-gauge perturbation. The SF in such a hybrid 
gauge would have the mode-sum formula 

oc 

^seif [hyb] = [^^d] tLA^-B"]~ SD" , (60) 

(=0 

where F^'' [rad] are the / modes of the full force in the radiation gauge, and the "counter 
term" SD" is the difference F^' [rad] — F^'' [hyb] , summed over / and evaluated at the 
particle. With a suitable choice of the hybrid gauge, the term 5Z)" can be calculated 
analytically, and Eq. ([60]) then prescribes a practical way of constructing the SF in a 
useful, regular gauge, out of the numerically calculated modes of the perturbation in 
the radiation gauge. 

Different variants of this idea were studied by several authors |1141 [5TJ 11441 [211 
1121) . but it has not been implemented in full so far. 

Generalized SF and gauge invariants. Another idea (set out in Ref. [21] and 
further developed in Ref. 77] ) involves the generalization of the SF notion through 
the introduction of a suitable averaging over angular directions. In some gauges which 
are not strictly regular in the aforementioned sense, the SF could still be defined in 
a directional sense. Such is the case in which the expression in square brackets in 
Eq. (j20p has a finite yet direction-dependent particle limit (upon transforming from 
the Lorenz gauge), and the resulting "directional" SF is bounded for any chosen 
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direction. (This seems to be the situation in the Regge- Wheeler gauge, but not in 
the radiation gauge — in the latter, the metric perturbation from a point particle 
develops a one-dimensional string-like singularity [H].) In such suitable 
averaging over angular directions introduces a well-defined notion of an "average" 
SF, which generalized the original MiSaTaQuWa SF (it represents a generalization 
since, obviously, the average SF coincides with the standard MiSaTaQuWa SF for all 
regular gauges, including Lorenz's). The notion of an average SF could be a useful 
one if it can be used in a simple way to construct gauge-invariant quantities which 
describe the radiative motion. This is yet to be demonstrated. 

A related method invokes the directional SF itself as an agent for constructing 
the desired gauge invariants. In this approach, one defines (for example) a "Regge- 
Wheeler" SF by taking a particular directional limit consistently throughout the 
calculation, and then using the value of this SF to construct the gauge invariants. This 
approach has been applied successfully, in combination with the mode-sum method, 
by Detweiler and others [S3l HI] ■ 

Radiation-gauge regular ization. Friedman and collaborators [96[ I95[ 167] 
proposed the following construction: Starting with the Lorenz-gauge S field, construct 
the associated gauge- invariant Newman-Penrose scalar -^q, and decompose it into 
spin- weighted spheroidal harmonics. Then obtain the harmonics of the retarded field 
ip'^Q^ by solving the Teukolsky equation with suitable boundary conditions, and (for 
each harmonic) define the R part through ijj^' = if/^^ — V'o • 1^ V'o known precisely, 
then ili^ is a vacuum solution of the Teukolsky equation. To this solution, then, apply 
Chrzanowski's reconstruction procedure to obtain a smooth radiation-gauge metric 
perturbation /i^^[rad], and use that to construct a "radiation gauge" SF (via, e.g., 
the mode-sum method). The relation between this definition of the radiation-gauge 
SF and the one obtained by applying the gauge transformation formula (|20p to the 
standard Lorenz-gauge MiSaTaQuWa SF is yet to be investigated. 

In reality, the S field is usually known only approximately, resulting in that -00 
retains some non-smoothness. How to apply Chrzanowski's reconstruction to non- 
smooth potentials is a matter of appreciable technical challenge. Also, the above 
procedure cannot account for the contribution to the SF from the two non-radiative 
modes, I — 0, 1, which then need to be treated separately, using other methods. 
Nonetheless, the technique offers a natural way of circumventing the gauge problem, 
and has much potential promise. 

Direct Lorenz-gauge implementation. In 2005, Barack and Lousto [19 succeeded 
in solving the full set of Lorenz-gauge perturbation equations in Schwarzschild, using 
numerical evolution in 1-l-lD. (An alternative formulation was later developed and 
implemented by Berndtson [33).) This development opened the door for a direct 
implementation of the mode-sum formula in the Lorenz gauge. It later facilitated the 
first calculations of the gravitational SF for bound orbits in Schwarzschild f25l ll431 [26]. 

In Sec. |6] we shall review this approach in some detail. Here we just point to 
a few of its advantages: (i) This direct approach obviously circumvents the gauge 
problem. The entire calculation is done within the Lorenz gauge, and the mode- 
sum formula can be implemented directly, in its original form, (ii) The Lorenz-gauge 
perturbation equations take a fully hyperbolic form, making them particularly suitable 
for numerical implementation in the time domain. Conveniently, the supplementary 
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gauge conditions (which take the form of eUiptic "constraint" equations) can be 
made to hold automatically, as we explain in Sec. [51 (iii) In this approach one 
solves directly for the metric perturbation components, without having to resort 
to complicated reconstruction procedures. This is an important advantage because 
metric reconstruction involves differentiation of the field variables, which inevitably 
results in loss of numerical accuracy, (iv) Working with the Lorenz-gauge perturbation 
components as field variables is also advantageous in that these behave more regularly 
near point particles than do Teukolsky's or Moncrief's variables. This has a simple 
manifestation, for example, within the 1+lD treatment in Schwarzschild: The 
individual multipolc modes of the Lorenz-gauge perturbation are always continuous at 
the particle, just like in the simple example of Sec. 14. II on the other hand, the multipole 
modes of Teukolsky's or Moncrief's variables are discontinuous at the particle, and so 
are, in general, the modes of the metric perturbation in the Regge- Wheeler gauge. 
Obviously, the continuity of the Lorenz-gauge modes makes them easier to deal with 
as numerical variables. 

5. 2. Numerical representation of the point particle 

Common to all numerical implementation methods is the basic preliminary task 
of solving the field equations (in whatever formulation) for the full (retarded) 
perturbation from a point particle in a specified orbit. This immediately brings about 
the practical issue of the numerical representation of the particle singularity. The 
particulars of the challenge depend on the methodological framework: In time-domain 
methods one faces the problem of dealing with the irregularity of the field variables 
near the worldline; in frequency-domain (spectral) treatments, such irregularity 
manifests itself in a problematic high-frequency behavior. We now survey some of 
the relevant methods. 

5.2.1. Particle representation in the time domain 

Extended-body representations. In the context of fully nonlinear Numerical 
Relativity, the problem of a binary black hole with a small mass ratio remains a difficult 
challenge, because of the large span of lengthscales intrinsic in this problem. (Current 
NR technology can handle mass ratios as small as 1 : 10 [76] — still nothing near the 
~ 1 : 10''-10^ ratios needed for LISA EMRI applications.) Bishop et al. [31] attempted 
a NR treatment in which the particle is modeled by a quasi-rigid widely-extended 
body whose 'center' follows a geodesic. Comparison with perturbation results did not 
show sufficient accuracy, and the method requires further development. 

An extended-body approach has also been implemented in perturbative studies. 
Khanna et al. [1011 [Wl [H] solved the Teukolsky equation in the time domain (i.e., 
in 2-l-lD, for each azimuthal m mode) with a source 'particle' represented by a 
narrow Gaussian distribution. This crude technique was much improved recently 
by Sundararajan et al. [152[ I153j using a "finite impulse representation" , whereby 
the source is modeled by a series of spikes whose relative magnitude is carefully 
controlled so as to assure that the source has integral properties similar to that of 
a delta function. Such methods were demonstrated to reproduce wave-zone solutions 
with great accuracy (indeed, that is what they are designed to do), but they are likely 
to remain less useful for computing the accurate local perturbation near the particle 
as required in SF calculations. 
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Extended-body representations suffer from the inevitable tradeoff between 
smoothness and locahzation: One can only smoothen the solution by making it less 
localized, and one can better localize it only by making it less smooth. In what follows 
we concentrate on methods in which the source particle is precisely localized on the 
orbit: The energy- momentum of the particle in represented by a delta- function source 
term [as in Eq. ([22]) ]. and the delta distribution is treated analytically within the 
numerical scheme, in an exact manner. 

Delta- function representation in 1+lD. In full 3-l-lD spacetime, the full 
(retarded) metric perturbation obviously diverges on the particle (at any given time) . 
The divergence is asymptotically Coulomb-like in the Lorcnz gauge (and can take 
a different form in other gauges). In spherically-symmetric spacetimes one can 
decompose the perturbation into tensor harmonics and solve a separated version 
of the field equations in 1+lD (time+radius) for each harmonic separately. In the 
particle problem this becomes beneficial not only thanks to the obvious dimensional 
reduction, but also because it mitigates the problematics introduced by the particle's 
singularity: The angular integration involved in constructing the individual I modes 
effectively "smears" the Coulomb-like singularity across the surface of a 2-sphere, and 
the resulting I modes arc finite even at the location of the particle. Furthermore, in the 
case of the metric perturbation in the Lorenz gauge, the individual I modes are also 
continuous at the particle [cf. Eq. ([50)1 in Seed]. The corresponding I modes of the 
Teukolsky or Moncrief gauge-invariant variables are generally not continuous at the 
particle, and generally neither are the modes of the metric perturbation in non-Lorenz 
gauges. 

The boundedness of the I modes is, of course, a crucial feature of the Z-mode 
regularization scheme, as we have already discussed. That same feature also greatly 
simplifies the 1-l-lD numerical treatment of the particle. Lousto and Price [ 106] 
formulated a general method for incorporating a delta-function source in a finite- 
difference treatment of the field equations in 1+lD. In this method, the finite-difference 
approximation at a numerical grid cell (in t-r space) traversed by the particle's 
worldline is obtained, essentially, by integrating the field equation "by hand" over 
the grid cell at the required accuracy. The original (ID) delta function present 
in the source term thereby integrates out to contribute a finite term at each time 
step. The original Lousto-Price scheme (formulated with a Ist-order global numerical 
convergence) was later improved by Martel and Poisson [108j (2nd-order convergence) 
and Lousto J103J (4th-order convergence) . This simple but powerful idea is at the core 
of many of the 1-t-lD finite-difference implementations presented in the last few years 
|107[ [80l [25] ■ including the work discussed in Sec. [6] below. 

Despite such advances, 1-l-lD particles remain numerically expensive to handle, 
because the nonsmoothness associated with them introduces a large scale variance 
in the solutions: The ^-mode field gradients grow sharply near the particle, and, 
moreover, become increasingly more difficult to resolve with larger I (recall the l- 
mode gradient is oc Z at large I). The mode-sum formula, recall, converges rather 
slowly (like ~ 1/0 1 ^-iid so requires one to compute a considerably large number of 
modes (typically ~ 20 with even a moderate accuracy goal). This proves to stretch 
the limit of what can be achieved today using finite differentiation on a fixed mesh. 

Several methods have been proposed to address this problem in the current 
context. Sopuerta and collaborators [1501 11491 [^ [Tfj explored the use of finite- 
element discretization. This technique is particularly powerful in dealing with multi- 
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scale problems, and, being quasi-spectral, it benefits from an exponential numerical 
convergence. So far it was applied successfully for generic orbits in Schwarzschild, and 
higher-dimensional implementations (for Kerr studies) are currently being considered. 
A related quasi-spectral scheme was recently suggested by Field et al. [BS]. Finally, 
Thornburg [157] very recently developed an adaptive mesh-refinement algorithm for 
Lousto-Price's finite differences scheme (with a global 4th-order convergence). This 
was successfully implemented for a scalar charge in a circular orbit in Schwarzschild, 
and generalizations are being considered. 

Puncture methods. In Kerr spacetime, one no longer benefits from a 1+lD 
separability. The Lorenz-gauge perturbation equations are only separable into 
azimuthal m-modes, each a function of t,r,9 in a 2-l-lD space. The m modes are 
not finite on the worldline, but rather they diverge there logarithmically (see the 
discussion in Sec. II. C of Ref. [17]). Since the 2+lD numerical solutions are truly 
divergent, a direct finite-difference treatment becomes problematic. However, since 
the singular behavior of the perturbation can be approximated analytically, a simple 
remedy to this problem suggests itself. 

The idea, which has recently been studied independently by several group 
[T71 11581 11051 1159j , is to utilize a new perturbation variable for the numerical time- 
evolution, which we shall call here the residual field. This field is constructed from 
the full (retarded) perturbation by subtracting a suitable function — the "puncture" 
field, given analytically — which approximates the singular part of the perturbation 
well enough that the residual field is (at least) bounded and differentiable at the 
particle. The perturbation equations are then recast with the residual field as their 
independent variable, and with a new source term (depending on the puncture field 
and its derivatives) which now extends off the worldline but contains no delta function. 
The equations are then solved for the residual field in the time domain, using (e.g.) 
standard finite differentiation. 

Several variants of this method have been studied and tested with scalar-field 
codes in 1-l-lD [158] and 2-l-lD [T7 ]fT05] . and also proposed for use in full 3-l-lD f95]. 
The various schemes differ primarily in the way they handle the puncture function 
far from the particle: Barack and Golbourn T7j introduce a puncture with a strictly 
compact support around the particle, Detweiler and Vega |158| truncate it with a 
smooth attenuation function, and in Lousto and Nakano I105j the puncture is not 
truncated at all. We will discuss the puncture method in more detail in Sec. [T] 

To obtain the necessary input for the SF mode-sum formula, the 2-l-lD (or 3-l-lD) 
numerical solutions need to be decomposed into I modes, in what then becomes a 
somewhat awkward procedure (we decompose the field into separate / modes just 
to add these modes all up again after regularization) . Fortunately, there is a more 
direct alternative: Barack et al. [28] showed how the SF can be constructed directly, 
in a simple way, from the 2 + ID m-modes of the residual field (assuming only that 
these mode are differentiable at the particle, which is achieved by designing a suitable 
puncture). This direct "m-mode regularization" scheme, too, will be described in Sec. 
[7| It is hoped that this technique could provide a natural framework for calculations 
in Kerr. It is yet to be applied in practice. 
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5.2.2. Particle representation in the frequency domain: the high-frequency problem 
and its resolution 

The I modes required as input for the SF mode-sum formula can also be obtained 
using a spectral treatment of the field equations. This has the obvious advantage 
that one then only deals with ordinary differential equations, although constructing 
the I modes involves the additional step of summing over sufficiently many frequency 
modes. (See [32] for a recent analysis of the relative computational efficiencies of 
frequency vs. time domain treatments.) As with the time-domain methods discussed 
above, the representation of the particle in the frequency domain too brings about 
technical complications, but these now take a different form. 

To illustrate the problem, consider the toy model of a scalar charge in 
Schwarzschild, allowing the particle to move on some bound (eccentric) geodesic of the 
background, with radial location given as a function of time by r = rp{t). Decompose 
the scalar field in spherical harmonics, and denote the multipolar mode functions by 
(I'lmityr). The time-domain modes are continuous along r — rp(t) for each l,m. 
However, the derivatives (f>im,r and (j)im,t will generally suffer a finite jump discontinuity 
across r = rp{t) (recall our elementary example in Sec. l4.ip . which reflects the presence 
of a source "shell" representing the I, m mode of the scalar charge. In particular, if 
the orbit is eccentric, the derivatives of (pim will generally be discontinuous functions 
of t at a fixed value of r along the orbit. 

Now imagine trying to reconstruct (pimit, r) (for some fixed r along the orbit) as 
a sum over its Fourier components: 



Since, for an eccentric orbit, (f>im{t,r) is only a C° function of t at the particle's 
worldline, it follows from standard Fourier theory [94j that the Fourier sum in Eq. 
(|6ip will only converge there like l/u. The actual situation is even worse, because 
for SF calculations we need not only but also its derivatives. Since (e.g.) (fum.r 
is a discontinuous function of t, we will inevitably face here the well known "Gibbs 
phenomenon" : the Fourier sum will fail to converge to the correct value at r ^ rp {t) . 
Of course, the problematic behavior of the Fourier sum is simply a consequence of our 
attempt to reconstruct a discontinuous function as a sum over smooth harmonics. 

From a practical point of view this would mean that (i) at the coincidence 
limit, r — > rp(t), the sum over lo modes would fail to yield the correct one-sided 
values of (l)im,r however many lu modes are included in the sum; and (ii) if we 
reconstruct (t>im,r at a point r = tq off the worldline, then the Fourier series should 
indeed formally converge — however, the number of lo modes required for achieving 
a prescribed precision would grow unboundedly as rp approaches rp{t), making it 
extremely difflcult to evaluate (f>im,r at the coincidence limit. 

This technical difflculty is rather generic, and will show also in calculations of 
the local gravitational field. The situation is no different in the Kerr case, because 
there too the mode-sum formula requires as input the sp/ier«ca/-harmonic modes of the 
perturbation field, and for each such mode the source is represented by a ^-distribution 
on a thin shell, which renders the field derivatives discontinuous across that shell. The 
problem becomes even more severe when considering gravitational perturbations via 
the Teukolsky formalism: Here, the l,m modes of the perturbation field (now the 
Newman-Penrose curvature scalar) are not even continuous at the particle's orbit — a 




(61) 
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consequence of the fact that the source term for Teukolsky's equation involves (second) 
derivatives of the energy- momentum tensor. Again, a naive attempt to construct these 
multipoles as a sum over their w modes will be hampered by the Gibbs phenomenon. 

A simple way around the problem was proposed recently in Ref. |29| . It was 
shown how the desired values of the field and its derivatives at the particle can 
be constructed from a sum over properly weighted homogeneous (source-free) radial 
functions Rimui{r), instead of the actual inhomogeneous solutions of the frequency- 
domain equation. The Fourier sum of such homogeneous radial functions, which are 
smooth everywhere, converges exponentially fast. The Fourier sum of the derivatives, 
which are also smooth, is likewise exponentially convergent. The validity of the 
method (and the exponential convergence) was demonstrated in Ref. [29] with an 
explicit numerical calculation in the scalar-field monopole case {I = 0). It was later 
implemented in a frequency-domain calculation of the monopole and dipole modes of 
the Lorenz gauge metric perturbation for eccentric orbits in Schwarzschild [511111431 
The same method should be applicable for any of the other problems mentioned above, 
including the calculation of EM and gravitational perturbations using Teukolsky's 
equation. 

The method of Ref. [29] (dubbed method of extended homogeneous solutions) 
completely circumvents the problem of slow convergence (or the lack thereof) in 
frequency domain calculations involving point sources. It makes the frequency-domain 
approach an attractive method-of-choice for some SF calculations. The method is now 
being implemented in first calculations of the scalar-field SF for Kerr orbits [162j . 

6. An Example: gravitational self-force in Schwarzschild via I+ID 
evolution in Lorenz gauge 

As an example of a fully worked out calculation of the SF, we review here the work by 
Barack and Sago on eccentric geodesic in Schwarzschild |25 1ll43j . This work represents 
a direct implementation of the mode-sum formula in its original form ([36]) . The 
decomposed Lorenz-gauge metric perturbation equations are integrated directly using 
numerical evolution in 1+lD. The numerical algorithm employs a straightforward 4th- 
order-convergent finite-difference scheme (a variant of the Lousto-Price method) on a 
fixed staggered numerical mesh based on characteristic coordinates. Below we briefly 
describe the perturbation formalism, discuss the numerical implementation in some 
more detail, and display some results. 

6.1. Lorenz-gauge perturbation formalism 

The linearized Einstein equations in the Lorenz gauge are given in Eq. (|22p . On the 
right-hand side of these equations is a distributional representation of the energy- 
momentum of a point particle of mass /i, which is moving on a timelike geodesic 
z^{t) of the background spacetime. Mathematically, the linear set ([22l) is a diagonal 
hyperbolic system, which admits a well posed initial-value formulation on a spacelike 
Cauchy hypersurface (see, e.g.. Theorem 10.1.2 of [161]). Furthermore, if the Lorenz 
gauge conditions of Eq. (fT4|l are satisfied on the initial Cauchy surface, then they 
are guaranteed to hold everywhere — assuming that the field equations are satisfied 
everywhere and that the energy-momentum satisfies V^'Tq/j = 0, as in the case of 
a geodesic particle. We remind that the gauge conditions do not fully specify 
the gauge: There is a residual gauge freedom within the family of Lorenz gauges. 
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hafj hap + V/3^ct + Vq^^, with any satisfying V"Va^^ = 0. It is easy to verify 
that both Eqs. (fTil) and remain invariant under such gauge transformations. 

Here we speciahze to eccentric geodesies around a Schwarzschild black hole, and 
employ a Schwarzschild coordinate system (t, r, 6, ip) in which the orbit is equatorial 
(z^ = 7r/2). Such orbits constitute a two-parameter family; we may characterize each 
orbit by the radial "turning points" rmin and rmax, or alternatively by the "semi-latus 

rectum" p = 2(rmax7'min)/(?'max + rmin) and "eccentricity" e = (rmax - rmin)/(rmax + 
rmin ) ■ 

Barack and Lousto jl9| decomposed the metric perturbation into tensor 
harmonics, in the form 

10 

hTp = 7 E E '^^^^""(-' ^) (62) 

and similarly for the source Tap. The harmonics Y^^^^"^ (0 , (p) (whose components are 
constructed from ordinary spherical harmonics and their first and second derivatives) 
form a complete orthogonal basis for 2nd-rank covariant tensors on a 2-sphere (see 
appendix A of [H]). The time-radial functions /i^*)'™ (i = 1, • • ■ , 10) form our basic 
set of perturbation fields, and serve as variables for the numerical evolution^ The 
tensor-harmonic decomposition decouples Eq. ((^ with respect to l,m, although not 
with respect to i: For each m, the variables Ji^*)'™ satisfy a coupled set of hyperbolic 
(in a 1-l-lD sense) scalar-like equations, which may be written in the form 

□ (2)j^«i™ + ^A/^WJJiO)'™ = S-W™ (i = 1, . . . , 10). (63) 

Here is the 1+lD scalar-field wave operator -I- V{r), where v and u 

are the standard Eddington-Finkelstein null coordinates, and V{r) — -1(1 — 
2M/r)[2M/r^ + l{l + l)/r'^] is an effective potential. The "coupling" terms 

AA^^^ji^h^^^''™ involve first derivatives of the /i'^^^'™'s at most (no second derivatives), 
so that, conveniently, the set decouples at its principal part. The decoupled 
source terms 5^'^'™ are each c>c 5[r — Vp^r)] (no derivatives of 6 function) and, as 
a result, the physical solutions /i^-')'™ are continuous even at the particle. Explicit 
expressions for the coupling terms and the source terms in Eq. (j63p can be found in 
Ref. [ig. 

In addition to the evolution equations ([551) . the functions /it*)'™ also satisfy four 
Ist-order elliptic equations, which arise from the separation of the gauge conditions 
(|14p into l,m modes. In the continuum initial-value problem, the solutions /i^*)'™ 
satisfy these "constraints" automatically if only they satisfy them on the initial 
Cauchy surface. This is more difficult to guarantee in a numerical finite-differences 
treatment, where (i) it is often impossible to prescribe exact initial data that satisfy 
the constraints, and (ii) discretization errors may amplify constraint violations during 
the numerical evolution. Inspired by a remedy proposed for a similar problem in 
the context of nonlinear Numerical Relativity [TS], Ref. [H] proposed the inclusion 
of "divergence dissipation" terms, oc V/iJ^^^, in the original set (P^ . so designed to 
guarantee that any violations of the Lorenz gauge conditions are efficiently damped 
during the evolution. These damping terms modify only the explicit form of the A4 
terms in Eq. ((63)) as shown in p9| . 

t To simplify the appearance of Eq. II62I I we have used here a normalization of h^*)''" which is slightly 
different from that of | 19| . 
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6.2. Numerical implementation 

The code developed by Barack and Sago [23 I143j solves the coupled set ([55)1 (with 
constraint dissipation terms incorporated in the M terms) via time evolution. The 
numerical domain, covering a portion of the external Schwarzschild geometry, is 
depicted in Fig. [3l The numerical grid is based on Eddington-Finkelstein null 
coordinates v, u, and initial data (the values of the 10 fields /i^')'™ for each I, m) 
are specified on two characteristic initial surfaces v =const and u=const. Equations 
(|63p are then discretized on this grid using a finite-difference algorithm which is 
globally 4th-order convergent. The numerical integrator solves for the various /i(*^'""s 
along consecutive u =const rays. A particularly convenient feature of this setup is 
that no boundary conditions need be specified (the characteristic grid has no causal 
boundaries). Moreover, one need not be at all concerned with the determination 
of faithful initial conditions: It is sufficient to set = on the initial surfaces 

and simply let the resulting spurious radiation (which emanate from the intersection 
of the particle's worldline with the initial surface) dissipate away over the evolution 
time. The early part of the evolution, which is typically dominated by such spurious 
radiation, is simply discarded. 



Figure 3. The numerical I+ID domain in the Barack-Sago code |25l 1143) : A 
staggered mesh based on characteristic (Eddington-Finkelstein) coordinates v, u. 
t and r* are, respectively, the Schwarzschild time and tortoise radial coordinates. 
The evolution proceeds from characteristic initial data on two null surfaces. 
Illustrated are a few sample geodesic orbits (radial plunge, circular, eccentric). 
[ Graphics reproduced from Ref. |12| .] 



The conservative modes I — and I — 1 (respectively the monopole and dipole) 
require a separate treatment, as they do not evolve stably using the above numerical 
scheme. (A naive attempt to evolve these modes leads to numerical instabilities which, 
so far, could not be cured.) Fortunately, the set simplifies considerably for these 
modes, and solutions can be obtained in a semi-analytic manner based on physical 
considerations. Detweiler and Poisson [55] worked out the / = 0, 1 Lorenz-gauge 
solutions for circular orbits, and their work is generalized to eccentric orbits in Ref. 
[50] . relying on the aforementioned method of extended homogeneous solutions. The 
calculation of [3D] yields the values of the fields /i^*)''" and their derivatives for Z = 0, 1. 
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To construct the I modes F^i_^ — the necessary input for the mode-sum formula 
— one first substitutes the expansion ([S^ in the expression for the retarded force 

[left-hand side expression in Eq. and then expands the result in spherical 

harmonics. The outcome is a formula for F^^[j_ given in terms of the various fields 

m ^j^j their derivatives (evaluated at the particle and summed over m). The 
number of Z'-modes /i^*)' contributing to each l-mode F^[^ depends on the off- 
worldline extension chosen for V"'^''. In [25 ) I143j we used a convenient extension in 
which [referring to Eq. ^ and Fig. [T] the metric g"^ takes its value at the field point 
X while the components u" (in Schwarzschild coordinates) take the fixed value they 
have at z. In this extension we find that the only contributing Z'-modes, for given 
I, are 1 — 3 < I' < 1 + 3. Explicit construction formulas for can be found in 

[^I143| . One then uses the numerical values of the fields and their derivatives, 

as calculated along the orbit using the above code, to construct F^[^ for sufhciently 
many I modes. The mode-sum formula p6l) then gives the SF. Figure 3] shows the 
final result for an eccentric geodesic with p — 7M and e = 0.2. 
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Figure 4. The gravitational SF [in units of (fi/M)^] along a Schwarzschild 
geodesic with semi-latus rectum p = 7M and eccentricity e = 0.2. The upper and 
lower lines show i^J^jf and F^"^^ , respectively. Integer values on the horizontal axis 
correspond to periapses (r = r^i^); note the slight retardation manifest in the 
magnitude of the radial component. The data for these plots were obtained using 
a direct integration of the metric perturbation equations in the Lorenz gauge, 
in conjunction with the mode-sum method, as described in Sec. |6] [Graphics 
reproduced from Ref. 1121 .] 

The calculation we have just described represents a milestone in the SF program: 
We now finally have a code able of tackling the generic EMRI-relevant SF problem in 
Schwarzschild. (Indeed, to address the fully generic physical problem would require a 
final generalization to the Kerr case!) Using this code, and similar codes developed 
independently by others [54j , we can now begin to explore the physical consequences of 
the gravitational SF — particularly those effects associated with its conservative piece. 
Work done so far includes (i) study of gauge invariant SF effects on circular orbits [54] ; 
(ii) comparison of SF calculations and results from FN theory in the weak field regime 
[Ml 1147] ; (iii) comparison of SF results from different calculation schemes and using 
different gauges [147j ; and (iv) a calculation of the shift in the location and frequency 



e = 0.2 
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of the Innermost Stable Circular Orbit (ISCO) due to the conservative piece of the SF 
pS] , Work is in progress to calculate SF-related precession effects for eccentric orbits. 
In Sec. [5] we will describe some of the recent work to explore the physical consequences 
of the gravitational SF. 

7. Towards self-force calculations in Kerr: the puncture method and 
m-mode regularization 

The time-domain Lorcnz-gauge treatment of Sec. [5] relies crucially on the separability 
of the field equations into (tensorial) spherical harmonics, which is no longer possible in 
Kerr. In the Kerr problem one can at best separate the metric perturbation equations 
into azimuthal m-modes, using the substitutior(|| 

oo 

E h':f,{t,r,9)e^^'^. (64) 

m— — oo 

Then one faces solving the coupled set for the 2+lD variables h™^. Insofar as 
vacuum perturbations are considered, this computational task is nowadays well within 
reach of even modest desktop computers. Indeed, over the past decade, 2+ ID 
numerical evolution has been a method of choice in many studies of Kerr perturbations 
[Ml ESI [Ml [911 mil ini [ing. The main challenge, rather, has to do with the 
inclusion of a 5-function source in the 2+lD numerical domain. This is problematic, 
because the 2+lD variables ft.™^ suffer a singularity at the location of the particle. 
The puncture method, which we have mentioned briefly in Sec. 15.2.11 overcomes 
this technical difficulty. In this section we review this method (as implemented 
by Barack and Golbourn [17 ) in some more detail. We also describe a new, ad 
hoc, regularization method — the m-mode regularization scheme — which enables a 
straightforward construction of the SF directly from the 2-l-lD numerical solutions, 
without resorting to a multipole decomposition. 

7.1. Puncture method in 2+lD 

In Sec. [2] we have spht the (trace-reversed) metric perturbation as Iv^^p — h^p + /i^^, 
with the gravitational SF then obtained from the smooth field h^^ as prescribed in 
Eq. Here we introduce a new splitting (defined below), 

so that 

h% = Kp + {TCr-hlp). (66) 

We denote the m modes of these various quantities, defined as in Eq. ([64]) . by /i^^™, 
'^c^'™' I16W splitting is defined (in a non- unique way) by introducing a 

puncture field /i^™'^, given analytically, which approximates h^p near the particle 
well enough that the m modes of the resulting residual field h'^^p are continuous and 

I A full separation of variables in Kerr is possible within Teukolsky's formalism, which, alas, brings 
about the metric reconstruction and gauge-related difficulties discussed in previous chapters. A full 
separation of the metric perturbation equations themselves, in Kerr, has not been formulated yet, to 
the best of our knowledge. 
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difFerentiable on the particle's worldline (and elsewhere). A particular such puncture 
is prescribed below [Eq. ([7S)l ]. The form of the puncture function away from the 
particle can be chosen as convenient — e.g., in such a way that it can be decomposed 
into m modes explicitly, in analytic form. 

The Lorenz-gauge perturbation equations (|22p are now written in the form 

V-^V./iLI + 2R'^c''phl'^ = Sl^^, (67) 

with 

^^'^l = -i6^r,;3 - v-<v,ICr - ^R^'phfTJ"'. (68) 

where — IGnTafs is the original (distributional) source term appearing on the right-hand 
side of Eq. ([22]) . The support of the source S^^^ now extends beyond the particle's 
worldline, but it contains no S function on the worldline itselfl^ The equations are 
separated into m modes, with the m-mode source given by 

Sl^;-'-' {t,r,0)^^ S^:i {t,r,0,^) e—'^ d^. (69) 

If the puncture is sufficiently simple, the source S'"^'™ can be evaluated in closed form 
(as in [H]). The m-mode field equations for the variables /i™|'™(t, r, 61), which are 
everywhere continuous and differentiable, can now be integrated numerically using 
straightforward finite differentiation on a 2-l-lD grid. 

To ease the imposition of boundary conditions for it is convenient to 

suppress the support of the puncture h^^'^ away from the particle, so that the 
physical boundary conditions for h^°p become practically identical to that of the 
full retarded field hap. In |17j this is achieved in a simple way by introducing an 
auxiliary "worldtube" around the particle's worldline (in the 2-f ID space): Inside this 
worldtube one solves the "punctured" m-mode equations for fi^^'"^, while outside the 
worldtube one uses the original, retarded-field m-modes /t™^'™ as evolution variables; 
the value of the evolution variables is adjusted on the boundary of the worldtube using 

7"ret,m Trcs,m , 7punc,m 

Two very similar variants of the puncture scheme have been developed and 
implemented independently by two groups jlTl I105j — both for the toy model of a 
scalar charge on a circular orbit around a Schwarzschild black hole (refraining from a 
spherical harmonics decomposition and instead working in 2-1- ID). The method is yet 
to be applied in Kerr and for gravitational perturbations. 

7.2. m-mode regularization 

The 2-l-lD numerical solutions obtained using the puncture method can be used to 
construct the input modes F^' for the mode-sum formula ((36|) , but this would require 
the additional step of a decomposition in spherical harmonics. It appears, however, 
that there is a simple way to construct the SF directly from the residual modes 
/i™^'™, without resorting to a multipole decomposition. Such "m-mode regularization" 
procedure was prescribed recently in Ref. [28 for the scalar, EM and gravitational SFs. 
We describe it here as applied to the gravitational SF. 

§ This can be shown by integrating S^'^J over a small 3- ball containing the particle (at a given time), 
and then inspecting the limit as the radius of the ball tends to zero 1171 . 
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In analogy with Eq. let us define the "force" fields 

F,ZM = = fiWf'^hl'^''^ (70) 

Then, recalling Eqs. and ([M]) . we may write the SF at a point z along the orbit 
as 

FZitiz) = hm [FZM + {F^nnM - Fs"(x))] . (71) 

Recall that the expression in square brackets {— fiW'^^'^h^^) is a smooth (analytic) 
function of x, and so the limit x ^ z is well defined. We proceed by prescribing a 
suitable puncture function ft-Jj""'^- 

As in previous sections, let us parameterize the particle's (bound, timelike) 
geodesic orbit by proper time r, and let z"(t) describe the orbit in Boyer-Lindquist 
coordinates. For an arbitrary spacetime point x" = (t, r, 9, ip) outside the black hole, 
let be the spatial hypersurface i=const containing x, let f{t) be the value of r at 
which Sj, intersects the particle's worldline, and denote z{t) = z[f{t)] and = u"(z). 
This procedure assigns to any point x outside the black hole a unique point z on the 
worldline, with four- velocity u°'. We denote the coordinate distance between a given 
point X and the point z associated with it by Sx"' = a;" — {t) , with Boyer-Lindquist 
components Sx" — (0, Sr, 59, Sip). 

Consider now the local form of the S field h^p{x, z), given in Eq. (I19|) . with x being 
an arbitrary point near the worldline, and z being the worldline point z associated with 
X. Recall in Eq. e{x; z) is the normal geodesic distance from x to the worldline, 
and u"(x; z) is the four- velocity parallelly propagated from z to x. The two quantities 
€^{x] z) and ■u"(x; z) are well defined and smooth functions of a; if a; is sufficiently close 
to the worldline [hence also close to z(t)]. For x ^ z (hence 5x 0) we have the 
asymptotic expansions 

e^(a;, z) = 5o + 5*1 + 0(fa*^), Ua{x,z)=Ua + 5ua + 0{5x'^), (72) 

where So and Si are, respectively, quadratic and cubic in 5x, and 6ua is linear in 6x. 
Explicitly, these expansion coeSicients read 

Sq = P^fjSx^Sx^ (73) 

Si =Pax^^Sx"Sx>'dx\ (74) 

dua = f^puxdx'^, (75) 

where Pajs = gap + gap = ga0{z) and f = r^^(z) (for a derivation of Si, see, 

for example, Appendix A of [22] )■ We now wish to regard Eqs. ([73l) -([75 |) as definitions 
of the quantities Sq, Si and 5ua, taken to be valid globally, for any x outside the black 
hole. [Of course, we keep in mind that the asymptotic expansions (17^ are only valid 
suflticiently close to the worldline.] We then define our puncture function, for any x, 
as follows: 

l^ap (5o + 5l)V2 • ^^^> 

The support of this function can be attenuated far from the particle in order to control 
its global properties, but such modifications will not affect our discussion here. 

It is not difficult to show [55] that, near the particle, the difference between the 
puncture ([75]) and the actual S field given in Eq. ([T5|) (with z taken to be z) has the 
form 

Cr(^) - = -^o'PapiSx) + const + OiSx'), (77) 
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where Eq = 5*0^^, and is a smooth function of the coordinate differences Jx", of 

homogeneous order (fe)*. [The function P^^^ depends on the function w^^ appearing 
Eq. HH), but here we wiU not need to know the exphcit form of these function.] It 
follows that the difference between the corresponding "force" fields has the local form 

F^^,,,{x) - Fiix) = e-Q 'P("5)(fe) + 0{Sx), (78) 

where P^^ is yet another smooth function, this time of homogeneous order (Sx)^. 
Notice that F^^^^ — Fg is bounded but generally discontinuous (direction-dependent) 
ai X ^ z. From Eq. (|7T|) it then follows that F^^ too is bounded but discontinuous 
at X ^ z [since the limit of the entire expression in square brackets in (j7ip is known 
to be finite and direction-independent]. 

We now arrive at the crucial step of our discussion. In Eq. (|7ip . for a given point 
z along the particle's worldline, we express the (analytic) function in square brackets 
as a formal sum over m-modes, in the form 

oo 

0(z)=lim ^ [F,",r(^) + (^^p"u™c(^)-^s"'"(^))]> (79) 

m — — OO 

where 

FZn^) C PrZsit. r, 9 ip'- z)e™('^-'^') dip', (80) 

and similarly for f^jj™ and F^'™" . Since the m-mode sum is formally a Fourier 
expansion, and since the Fourier expansion of an analytic function converges uniformly, 
we may replace the order of limit and summation in Eq. ([79]): 



FU{z)= E lim[F,",r(^) + (^^p"u"c(^)-^s"'"(^))]- (81) 



X — yz 

m— — oo 



From Eq. fTSl) . omitting the terms 0{5x) as they cannot possibly affect the final value 
of the SF in Eq. jTlI), we have 



hm {F^-Z-Ft"') = lim T e^' P^.-^e-^^^^' d^' , 



(82) 



where the integrand is evaluated a,t (p = (p' and where we have used the fact that 
X — > z implies also x —^ z. Crucially, one finds |28] that the last integral vanishes at 
the limit x ^ z, for any m, regardless of the explicit form of ^'(5)- Hence, Eq. ([5T|) 
reduces to 

00 

^^"elf(^)= E ^rcr(^)- (83) 

m— — 00 

Here the substitution a; = z is allowed since the limit x —i- z is known to be well 
defined (and, in particular, direction-independent). Phrased differently, the modes 
F^'J" corresponding to our puncture (ffS)) are continuous at the particle, as desired. 

It might perhaps seem suspicious that, in the above analysis, the sum over m- 
modes F^^^^ — Fg'™ (which are all zero at x = z) fails to recover the original function 
■^punc ~ ^s' (which is discontinuous at x = z and hence indefinite there). This, 
however, should not come as a surprise. Recall that the formal Fourier sum at a 
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step discontinuity (where the function itself is indefinite) is in fact convergent: it 
yields the two-side average value of the function at the discontinuity. Technically, this 
peculiarity of the formal Fourier expansion is due simply to the non-interchangeability 
of the sum and limit at the point of discontinuity. 

Relatedly, we should emphasize that the full (4D) residual field, i^j-esj does not 
yield the correct SF upon taking the limit x z: this field does not even have a well 
defined limit a; ^ z, as we argued in the discussion below Eq. ((78|l . However, the 
sum over the formal Fourier modes of F^^, which indeed fails to reproduce F^^ at the 
particle, does turn out to give the correct SF, as we have established in the discussion 
leading to Eq. 

It is possible to design an improved, higher-order-accurate puncture h^^'^, so 
that the SF is indeed simply given by F^^{z). Such an improvement is necessary 
for 3-l-lD implementations of the puncture scheme, and is beneficial also in a 2-l-lD 
framework as it enhances the differentiability of the residual variable and improves the 
convergence of the m-mode sum in Eq. (j83p . However, such improvements entail an 
explicit calculation of higher-order terms in the S field [including the term w^p and 
possibly higher-order terms in Eq. (fT9|) ]. A higher-order puncture, suitable for 3+ ID 
implementation, was devised recently by Vega and collaborators [158|,I159] . 

Our result ([83|) can be written explicitly in terms of the m modes of the residual 
field h'f^^, and easily so if we use the fixed off-worldline extension of V"'^'*' (the choice 
of extension may affect the value of the individual m modes, but not the eventual 
value of the SF). We then have 



where, of course, the right-hand side is evaluated at the particle. This formula 
prescribes a straightforward method for constructing the SF in a 2-1- ID framework. 
In the puncture scheme we effectively "regularize" the field equations themselves (not 
their solutions, as in the standard Z-mode regularization method), by writing them 
in the form ()67|) with a sufficiently accurate puncture function — like the one we give 
in ([7^ . Once the m- modes of the residual perturbation have been calculated, the 
SF is constructed directly from these modes, via Eq. ([M]) . with no need for further 
regularization. This "m-mode regularization scheme" , yet to be applied in actual 
calculations of the SF, offers a simple and efficient framework for such calculations in 
Kerr spacetime. 

8. Physical effects of the self force 

With the development of the first gravitational SF codes in the past few years, it 
became possible to start exploring quantitatively the physical consequences of the SF. 
While the ultimate goal of the SF program remains the description of the long-term 
orbital evolution, knowledge of the SF along geodesic orbits already gives access to 
some interesting physics associated with the 0(/i) dynamics. Among the questions 
that one can address: How does the finite mass of the particle affect the rates of 
periapsis precession and orbital plane precession? How does it modify the location and 
frequency of the ISCO? Are there other, less familiar conservative SF effects that could 
manifest themselves in a characteristic way in the emitted gravitational waveforms? 



oo 




(84) 



m— — OO 
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Can transient resonances between the radial and polar motion (for eccentric and 
inclined orbits in Kerr) have important observational implications [55] ? 

Identifying and quantifying concrete SF effects that are gauge invariant and have a 
clear physical interpretation is also vital if one intends to cross- validate SF calculations 
carried out in different gauges |147| , and in making a connection with results from PN 
theory [54|- Quantitative information about concrete SF effects [such as the 0{fi) 
correction to the strong-field periapsis precession rate as a function of the geodesic 
parameters] can be incorporated "by hand" into approximate (PN) models of EMRI 
orbital evolution [20] ■ In this way, the study of SF effects for geodesic orbits can 
inform the improvement of EMRI models even before a reliable and workable scheme 
for the long-term evolution is at hand. 

In this penultimate section we describe some of the recent initial investigations 
into the physical consequences of the gravitational SF. 

8.1. Conservative and dissipative pieces of the SF 

In discussing the physical consequences of the SF it is very useful to distinguish 
between "conservative" and "dissipative" effects. These are normally defined as effects 
arising, correspondingly, from the "time symmetric" and "time antisymmetric" pieces 
of the SF [SJ. To formulate this more precisely, let us re- write Eq. ([M]) as 

O(rct) (^) = lim [FZt {x)-F^{x)], (85) 

where the label 'ret' is to remind us that the physical SF is derived from the physical, 
retarded metric perturbation [recall F^"^ — fJ-'^x^'^f^'p^]- Then introduce the force 

O(adv) (^) ^ lim [F:^. (x) - F^ (x)] , (86) 

where F^^^v = t^-^Z'^'^hf:;' (namely, ^;eif(adv) is obtained from F^^^^^^^^^^ by replacing 
the retarded perturbation with the advanced one). The conservative and dissipative 
pieces of the SF are then defined through 

-^cons = 2 (^SGlf(rct) + ^SGlf(adv)) ' -^<?iss = 2 ('^sclt(rct) ~ ^SGlf(adv)) • i^'^) 

The physical SF is the sum of the two pieces: F^^,, = F^^,,^^^^^ = F<?„„, + F^^^. 

The dissipative force, F^-^^^ is responsible (in particular) for the long-term secular 
drift in the value of the intrinsic "constants" of motion — the energy azimuthal 
angular momentum and Carter constant Qc associated with the geodesic motion 
in Kerr. -Fc'onsj on the other hand, has no such long-term influence on the evolution 
of these orbital parameters. Here "long-term" refers to a period of time much longer 
than the longest orbital periodHI It should be noted, however, that both F^^^^ and 
■^cons affect the momentary values of the intrinsic parameters E,Lz,Qc (in a gauge- 
dependent way); in the case of fc'ons this effect "averages out" over many orbital 
periods, whereas for F^^^^ it accumulates. It should also be noted (what is usually 
less well recognized) that F^^^^, too, gives rise to secular, long-term effects — those 
associated with the evolution of the phase-type orbital elements [ 1351 1136[ 1137] . 



t Recall bound geodesies in Kerr are 3-periodie, with three characteristic frequencies corresponding 
to the azimuthal, radial, and longitudinal motion. 



L.Barack 



42 



How does one go about separating the full SF into its conservative and dissipative 
pieces in practice? In a frequency-domain framework this poses no difficulty: In 
addition to calculating the physical SF Fseif(rct)j one also calculates the advanced 
metric perturbation modes by suitably inverting the boundary conditions in the 
relevant ordinary differential equations, and use these modes to construct the quantity 
-^soif(adv) b-g-' using the mode-sum formula ([361), replacing F^l^ ^rdv±]- The 
conservative and dissipative pieces are then obtained using Eq. (j87p . This procedure 
was implemented, for example, in the analysis of Ref . [5D] . The situation is slightly 
different in a time-domain framework, as one does not normally control the boundary 
conditions during the time evolution (cf. Sec. [B])- Haas [81 proposed to obtain the 
necessary advanced SF data by reversing the time-direction of the evolution: For a 
given orbit, one solves the relevant time-domain field equations once evolving forward 
in time as usual to obtain -Fscif(rct)i ^^id then again evolving backward in time (starting 
with initial Cauchy data specified in the future) to obtain i^sGif(adv) • 

There is a more computationally economic method for constructing ^'(.'ons/diss 
the time domain, at least in the case of orbits that are either equatorial or circular. 
This method, first proposed and implemented in Ref. [26], makes use of the symmetries 
of Kerr geodesies, first noted in the current context by Mino |110j . It can be shown 
(see, e.g.. Sec. II. G of [87j) that, at any given point in Kerr spacetime, the following 
relation applies: 

O(adv)("t''«r,W0,Uv) = Ho')Kelfiret)iut,-Ur,-Ug,U^), (88) 

where 

= (-1,1, 1,-1) (89) 

in Boyer-Lindquist coordinates (no summation over a on the right-hand side). Here 
we are treating the SF as a function of the four- velocity at the given spacetime point. 
Now consider a bound geodesic orbit in Kerr, which is either circular (and possibly 
inclined) or equatorial (and possibly eccentric) . We parameterize this geodesic by the 
proper time r, and, if the orbit is eccentric, take r = at one of the radial "turning 
points" (i.e., where dr/dr — 0) without loss of generality. Now examine the retarded 
and advanced SFs ^^scif(rct/adv) (''') S'lo'^g this geodesic. It is clear from symmetry that 

^solf(rot/adv)K''«r,Me,M>p;-T) = F°if(„t/^dv) ("* ' ""r ' ""e , ; t) . (90) 

From Eq. dH]) it then follows that 

'^solf(rot/adv) (''') = '^(a) -^SGlf(adv/rct) ("''') ' (^l) 

and Eq. ^ gives 

^cons(T) = I (0(rct)(^) + ^ (a) O(rct) (-^)) : (92) 
^d"ss(^) = I (0(ret)(^) - ^(a)0(rct)(-^)) • (93) 

Since in time-domain calculations the physical (retarded) SF is computed along an 
entire radial period (at least), Eqs. ([5^ and can be used to obtain the conservative 
and dissipative components of this force anywhere along the orbit without resorting 
to a calculation of the advanced perturbation. 

The above trick can, of course, be implemented for any geodesic orbit in 
Schwarzschild (as in this case the orbit can always be taken as "equatorial"). However, 
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it cannot be applied for orbits in Kerr that are both eccentric and inclined, because 
for such orbits the symmetry relation (^0]) does not hold in general. Note that, for 
circular orbits, Eqs. and ([55)1 immediately yield F^^s = -^cons = as well as 
-^diss — ^diss — 0, so in this case the t and ip components of the SF are purely 
dissipative, while the r and 6 components are purely conservative. 

Finally, it is useful to have at hand separate mode-sum formulas d la (|36p for 
the conservative and dissipative pieces. To obtain such formulas, first write a mode- 
sum expression for ^;eif(adv) by replacing F^"J^± F^j^^ in Eq. Here F^^^vi 
are derived from the modes of the advanced metric perturbation in just the same 
way as F^eti derived from the modes of the retarded perturbation. (The same 
regularization parameters apply to both retarded and advanced forces, because the 
corresponding metric perturbations share the same singular piece.) Then substitute 
the mode-sum expressions for Fg"if(rgt) and Fgoif(adv) f^T]) . This gives 

oc 

^cons = E (^co„.± T LA^ - B^) , (94) 

oo 

^diss = E^dt^ (95) 

(=0 

where we have used C" = L»" = and introduced i^eonsi = (-P'rcli + ^a"dv±)/2 and 
-^dissi = (-^roti ~ ^iSiv±)/2- Noticc that the dissipative component of the SF requires 
no regularization. In fact, one can show that the mode sum in Eq. ()95|) converges 
exponentially fast. For that reason, the computation of the dissipative piece of the 
SF via the mode-sum method is technically much less challenging than that of the 
conservative piece. 



8. 2. Dissipation of energy and angular momentum 

Perhaps the most familiar aspect of the self-interaction physics in the binary context 
is the long-term radiative decay of the orbit. In the language of the SF, we say 
that the work done on the particle by the dissipative component of the SF converts 
orbital energy (and angular momentum) into gravitational-wave energy (and angular 
momentum). The relation between the SF and dissipation is immediately evident 
from the equation of motion (HI]), whose (Boyer-Lindquist) t and (p components read, 
respectively, 

^^ = Fr^ /^^ = ^r- (96) 

In the absence of SF (that is, in the test-particle limit), Eqs. ((96)) tell us that ut and 
Uip are constants of the motion, and we interpret these as the orbital energy, E = —ut, 
and azimuthal angular momentum, Lz = — both per unit charge Equations 
(|96p also tell us how E and Lz change under the effect of the SF. This effect is not 
entirely dissipative: The SF includes periodic components which do not lead to a net 
long-term change in the values of E and Lz- The non-periodic component of Ft and 
F^ does lead to a secular drift in the values of E and L, interpreted as dissipation. 

This non-periodic component is entirely contained in the dissipative piece of F^ 
and F^^ as expected. This is immediately clear from Eq. ([^^ in the case of equatorial 

X This interpretation is motivated by the observation that — ut and Uip are, in fact, projections of the 
four-velocity onto the Killing vectors associated with the invariance of the Kerr background under 
(respectively) time translations and spatial rotations. 
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or circular orbits in Kerr (and for all orbits in Schwarzschild): The anti-symmetry 
of F™"''(t) and F™"''(t) with respect to r — > — t means that these pieces of the SF 
vanish upon integration over one complete (radial or longitudinal) period, and thus 
produce no net change in the values of E and over that period (this, of course, 
assumes that the orbit is very nearly geodesic over one period). The same applies 
generically for all orbits in Kerr, if the orbital integration is taken over sufRciently 
many periods (formally, infinitely many). These statements are encapsulated in the 
relations 

^,{E)^~{Ft/u')=-{Ff'^yu'), fi(L,) = {Fju')^{F^'^yu'), (97) 

where an overdot represents d/dt (hence the factor = dr/dt on the right-hand 
side), and (•) indicates an average over sufficiently long time. One always finds 
^^diss^ > and (F^*^") < (see, e.g., [122]), so that {E) < and (L,) < 0, and 
orbital energy and angular momentum are lost over time as one expects. The time- 
average quantities in Eq. (j97p are expected to be independent of the gauge even though 
the momentary values of F^^^^ and F^^^^ are themselves gauge dependent. 

Implicit in the above discussion is the assumption that the orbit is indeed 
"available" for the necessary time-averaging, i.e., that E and Lz evolve very slowly over 
the required averaging time — otherwise the averaging procedure would be meaningless. 
Such "adiabaticity" requirement (formulated more accurately, e.g., in [91j ) is believed 
to hold well in LISA-relevant EMRIs throughout the inspiral and until very close to 
the innermost stable orbit — see, e.g., [91j for a more quantitative discussion of this 
point. 

The loss of orbital energy and angular momentum is precisely balanced by the 
flux of the corresponding quantities in the gravitational waves radiated to infinity and 
down into the black hole. Using Eq. ([57)1 we can express this balance in terms of the 
local SF: 

-f,-\F^'^yu') = (i?w>eh + (i?^)oo. (98) 

The quantities on the right-hand side are the time-averaged asymptotic fluxes down 
the event horizon ('eh') and out to infinity ('oo'), with the convention that the infinity 
fluxes are positive. This convention also fixes the sign of the horizon fluxes for a given 
orbit, but it should be noted that the latter can be either positive or — for certain 
orbits in Kerr — negative. Negative horizon fluxes reflect a superradiance behavior, 
whereby energy and angular momentum are, in effect, transferred from the ergosphere 
of the Kerr hole to the orbit | 1561 [9T| [75] . The formal proof of the balance equations 
(|98p involves the application of the Gauss theorem in the 3-volume bounded between 
a 2-sphere at r — > oo and another one just outside the horizon; how this is done is 
demonstrated, for example, in appendix D of Ref. [54j (in the special case of circular 
orbits in Schwarzschild). 

The asymptotic fluxes (£''^^)oo and (i^^)oo are calculated across a large 2- 
sphere r oo residing in the "wave zone" , where the background's radius of curvature 
is much larger than the characteristic wavelength of the emitted gravitational waves. 
In this case Isaacson's effective energy-momentum tensor formulation [93| applies, and 
one can use it to calculate the fluxes at infinity. This is a standard calculation in black 
hole perturbation theory, and since the early 1970's it has been performed many times 
using frequency-domain methods (e.g, [1561 11301 [OTl [75] ) and more recently within 
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a time-domain framework |107l 1191 I25| . There is also a standard prescription for 
calculating the horizon fluxes {E'^'^)oc and (ij^)oo; it is due to Teukolsky and Press 
|156| (based on the horizon perturbation formalism of Hawking and Hartle [55) and 
is formulated in the frequency domain in terms of curvature scalars. A time-domain 
formulation of the horizon fluxes was more recently developed by Poisson [134j . 

The balance equations ((98|) . with ([QT]) . allow us to infer the leading-order long- 
term evolution of E and Lz without resorting to an explicit calculation of the local 
SF (assuming adiabaticity). To get a full description of the time- average orbital decay 
one must, in general, also be able to calculate the evolution of the Carter constant 
Qc- It is not difficult to write down a relation analogous to ([97]) which expresses 
(Qc) directly in terms of the local SF [122j . but there is no known analogue to (jOS]) . 
relating (Qc) to the asymptotic flux of radiation. However, thanks to a breakthrough 
idea by Mino |110| llllj and follow up work by Sago and collaborators |146| [73] , there 
is now known a practical formula for {Qc) which does not require knowledge of the 
SF (it does require knowledge of the local advanced and retarded modes of the metric 
perturbation) . 

In conclusion, information on the time-averaged evolution of all three intrinsic 
"constants of motion", E, Lz and Qc, is directly accessible from the dissipative piece 
of the SF, but there are alternative methods (likely more computationally efhcient) 
for accessing this information without resorting to the SF and to Eq. (|98p . However, 
from the perspective of the SF program development, the balance equations ([98]) are 
still a very useful tool for validating SF codes. Two independent tests are possible: 
First, since a SF code always involves a calculation of the metric perturbation itself, 
one can use this perturbation to construct the asymptotic fluxes of energy and angular 
momentum, and then use Eqs. (|55]) to test the self-consistency of the code. Second, 
one can test the computed values of the dissipative SF against asymptotic flux data 
available in the literature (e.g., |107[ [75]). Such a two-fold test is an important, 
reassuring check on the validity of a SF code. Explicit calculations demonstrating the 
balance relations ([98]) have been carried out so far for circular 25] and eccentric [143] 
orbits in Schwarzschild. 

8.3. Conservative effects on circular orbits in Schwarzschild 

Detweiler was first to explore quantitatively the conservative effects of the SF — in 
the simple example of circular geodesies outside a Schwarzschild black hole ^54]. He 
pointed out two gauge-invariant quantities that carry non-trivial information about 
the conservative SF dynamics in this setting: the orbital frequency (with respect to 
time t), 

n-u'^/u^ (99) 
and the contravariant t component of the four-velocity, denoted 

U = u\ (100) 

Here — dx°^ jdr denotes the value of the perturbed four- velocity along the circular 
geodesic orbit, including SF corrections. We have m'' = along the orbit, and we also 
take (without loss of generality) = 0. The said gauge invariance of and U is 
restricted to transformations for which the gauge displacement vector ^" respects the 
helical symmetry of the perturbed spacetime system, i.e., it satisfies 



(101) 
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through 0{fi), at least along the particle's worldline. It is easy to see that all 
contravariant components are invariant through 0{ii) under such transformations 
|147| . The two nontrivial components u* and give two independent gauge invariants, 
which we can re-conibine to form Detweiler's gauge-invariant variables Q and U . 

While the physical significance of the frequency is clear, that of U is slightly 
less obvious. Detweiler |54] discusses two physical interpretations of U. First, it is a 
measure of the gravitational red-shift experienced by photons emitted by the orbiting 
particle and observed at a distance. Second, U is intimately related to the helical 
Killing vector of the perturbed spacetime [with its gauge invariance being simply the 
statement that this Killing vector is invariant under gauge transformations satisfying 
Eq. ((99|) ]. Unfortunately, as we explain below, the conservative SF corrections to Q 
and U are not expected to have any short-term observable imprint on the gravitational 
waveforms emitted from the circular orbits. 

Explicit expressions for fl and U, including SF terms, are easily obtained from 
the r component of the equation of motion (HJ), setting = and du^ /dr = 0. One 
finds [143 

ro(ro - 3M) 



2nM 



(102) 



U = Uo(^l- ^Frj , (103) 

through 0(fi), where rg is the orbital radius (Schwarzschild r coordinate), and 
fio = (M/r^)^/^ and Uq = {1 ~ 3Af/ro)^^/^ are the geodesic (unperturbed) values 
of and U, respectively. Recall that the radial component of the SF, F^, is purely 
conservative in the case of a circular orbit; hence the SF terms in Eqs. (|102p and 
(|103p represent conservative corrections to ft and U. Of course, the two quantities 
would also evolve dissipatively (this effect is described by the t and tp components of 
the equation of motion), but here we ignore the dissipative piece of the SF in order 
to study the effect of its conservative piece in isolation. Equations (|102|1 and (|103p 
tell us that the effect of the conservative SF is to "shift" the values of fl and U as 
compared with their non-perturbed geodesic values (for a given orbital radius tq). 
The conservative SF shift in 17, AH, = fl — Hq, was first calculated (numerically, as a 
function of tq) in Ref. [25], along with the conservative shift in the orbital energy E 
(which is gauge dependent). 

It should be understood that, despite the formal gauge invariance of fl, the 
relation AO(ro) is, in fact, gauge dependent, because the radius Tq itself is gauge 
dependent (the calculation in |25| was carried out in the Lorenz gauge). To illustrate 
this point, suppose that we carry out two independent calculations of the SF in two 
different gauges, and we wish to test our results by comparing the values of using 
Eq. (|102p . The problem we would face is that a given value of ro would, in general, 
correspond to two physically distinct orbits, due to the gauge ambiguity at 0{fj,); 
there is no way of relating the coordinate values of the two physical radii in the two 
gauges without further information about the gauge (in the form of the local metric 
perturbation, for example). The conservative SF shift in U, AU{ro) = U — Uq, is 
similarly gauge dependent, and thus it too cannot be utilized usefully as a gauge- 
invariant measure of the SF effect. 

One might hope to get around this problem by expressing one of our gauge 
invariants in terms of the other. To this end, it is convenient to introduce the gauge 



L.Barack 



47 



invariant radius, 

i?= (M/^72)'/^ (104) 
and then express U in terms of R. However, one then simply finds 

U{R) = {1-3M/Ry^^^ + 0{n'^), (105) 

which contains no exphcit information about the SF. This is an obvious result: An 
0(/i) term on the right-hand side of Eq. (|105p could only involve R and Fr, and since 
R is gauge invariant while is not, the occurrence of such a term would necessarily 
violate the gauge invariance of U. We therefore have to maintain our conclusion that 
U and Q (or R), each on its own or even both combined, do not contain gauge- 
invariant information about the SF effect, unless, in addition, we have access to local 
gauge information. For that reason, also, we cannot expect to be able to identify the 
effect of the conservative SF in a detected gravitational wave from a circular EMRI 
merely by measuring the momentary values of and/or U (even if there was a way 
of extracting U from the waveform, which is unhkely)I3 

Fortunately, information about the local metric perturbation, not available to 
the gravitational-wave astronomer, is available to the SF theorist running a SF code. 
This allows the theorist to utilize and U usefully in quantifying the gauge-invariant 
content of the conservative SF, as first shown by Detweiler [S31 [SI] and further 
demonstrated in [147j . In the following we outline the method of Ref . [147] . 

Consider again a circular orbit of radius tq (in Schwarzschild) , but now interpret 
this, a la Detweiler and Whiting, as a geodesic of the effective perturbed spacetime 
with metric gajs + h^^, where h^^ is the i?-part of the physical perturbation associated 
with the particle. Let f be proper time along this geodesicU To a given physical event 
along the orbit there thus correspond two proper time values, r and f, associated with 
the two interpretations. (The two values will differ, in general, even if we calibrate 
T and f to agree with one another at some initial moment.) It is easily shown |147j 
that, through O(^), 

dr 1 

— = 1 + H, with H = -h^gu"u>^, (106) 
dr 2 

where the perturbation is, of course, evaluated at the worldline point in question. [The 
four-velocity in H can be defined interchangeably with respect to either r or f: as 
already 0{fj,), the difference would only affect terms of 0{ii^) which are anyway 
neglected in our discussion.] Next, re-define the gauge invariants U and Vl in terms of 
the four- velocity = dx"^ /dr: 

fi = u^/w* = u^/u* = f], C/ = u* = t/(l + iJ), (107) 

where we have used Eq. (jl06p . Finally, express U in terms of the gauge invariant 
radius R^ R, and construct the 0{y) difference MJ{R) = U{R) - (1 - iM/R)-^/"^. 
The final result is [HT] 

AU{R) = (1 - MI/R)-^'^H. (108) 

f It is less clear how the conservative piece of the SF might affect the long-term evolution of the 
circular orbit under the full SF (dissipation included) , and what influence it may have on the emitted 
gravitational waveform — these questions await investigation. 

I Note our notation here for t and f is reversed compared to that of Ref. |147) . 
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Evidently, the metric perturbation combination H is gauge-invariant, as can be verified 
with an explicit calculation. 

Equation (jlOSp provides a nontrivial gauge-invariant relation which explicitly 
involves the i?-part of the local metric perturbation (although it does not involve 
directly the SF). Two theorists working in two different gauges should be able to 
agree on the value of AU{R), and such an agreement would constitute a nontrivial 
test of the calculation of h^^ — and, to an extent, of the SF too. 

Two such tests, based on Eq. (|108p . were carried out so far. In Ref. jl47j the 
results of Barack-Sago's Lorenz-gauge SF code [25j were compared with those of 
Detweiler's Regge- Wheeler SF code [Sj. The two codes were shown to agree on the 
values of AU{R) to within the computational error (of merely ~ 10^^ in fractional 
terms). In Ref. |54], Detweiler worked out a 2PN expression for AU{R) and compared 
it with the numerical SF data — showing an astonishing agreement down to radii as 
small as i? ~ 8M. Later, a 3PN expression derived by Le Tiec and collaborators |100j 
showed an even closer agreement. The 3PN expression approximates the "exact" SF 
value of At7 to within mere ~ 1% at i? = 12M and ~ 5% at i? = 7M. 

8.1 IS CO shift 

To identify conservative SF effects with truly "observable" consequences (by which we 
mean ones that can be measured in the gravitational waveform at least in principle) , 
we must move away from the simplicity of circular orbits. We need not move very 
far, though. There is a simple effect with a clear physical significance, which is 
manifest already in the dynamics of orbits with infinitesimally small eccentricity: the 
conservative SF-induced shift in the value of the orbital frequency 51 at the ISCO. 
This frequency shift, which is gauge invariant [under transformations satisfying Eq. 
(|10ip ]. was calculated recently in Ref. and we shall describe this analysis briefly 
below. 

Strictly speaking, the ISCO is defined in a precise way only at the test particle 
limit, /z — > 0. When /i is finite, dissipation "smears" the transition from inspiral 
to plunge and it is no longer precisely localizable. Ori and Thorne showed |124j (for 
circular orbits in Kerr, a result later generalized to other orbits by O'Shaughnessy [125] 
and Sundararajan |151| ) that the "width" of the radiative transition regime (measured, 
e.g., in terms of the frequency bandwidth of the emitted gravitational wave) scales 
as a low power of the mass ratio: ^ {^/MY^^ . The same scaling was discovered 
independently by Buonanno and Damour for a binary of nonrotating black holes with 
an arbitrary mass ratio [35] ■ Ori and Thome's analysis ignores the conservative effect 
of the SF, but the latter is expected to modify the location of the ISCO only by 
an amount ^ (/i/M), which for very small mass ratios we expect to be negligible 
compared with the width of the radiative transition. Still, there is a strong motivation 
to study the conservative effect of the SF on the ISCO. First, we would like to confirm 
the above expectation and quantify it better. Second, the conservative 0(/i) shift 
of the ISCO frequency (dissipation ignored) is a precisely specifiable gauge-invariant 
quantity, and as such it can serve as a convenient "anchor" point for comparison 
between different calculation schemes. In particular, one can envisage it being used 
as a strong-field benchmark for calibration of FN calculations. Indeed, the value of 
the conservative ISCO frequency shift, for mass ratios not necessarily small, has been 
utilized extensively in the past by PN theorists in assessing the performance of various 
PN schemes (see, e.g., |51[|35lf37] ). The calculation in Ref. [26 now provides this value 
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precisely, at 0{fj,), in the Schwarzschild case. 

To outline the calculation, we first remind how the occurrence of an ISCO is 
explained from the point of view of an effective radial potential. Ignoring the SF for 
the moment, the radial component of the equation of motion for a timelike geodesic 
in Schwarzschild takes the familiar form 

/i dV^g{r,L,) _ 
/^^ = ~2 0^ = -^<=ff' (10^) 

where = (1 — 2Af/r)(l + L^/r^) is an effective potential for the radial motion, and 
throughout our current discussion an overdot will denote d/dr (not d/dt as in Subsec. 
18. 2p . The quantity J-cg{r, Lz) can be interpreted as an effective radial force acting 
on the geodesic test particle. A stable circular orbit is associated with the (single) 
minimum of Vcs (for given L^), which occurs only for > \/l2M. The radius of 
the circular orbit, r — tqILz), decreases monotonically with decreasing L^, and at the 
limiting value of Lz = vl2M we have ro = 6M — the innermost stable circular orbit. 

To better understand the dynamical significance of the ISCO it is instructive 
to examine the behavior of a circular orbit under a small-eccentricity perturbation. 
Writing r(T) — + e6r(T) and considering the linear variation of Eq. (|109p with 
respect to the small eccentricity e, one readily obtains [26] 

^r(r) = -w^(5r(r), (110) 

with 

^ M(.o-6M) 
° r3(ro-3M) ^ ^ 

Thus the radius of the perturbed orbit is a linear oscillator with frequency ujq [the 
subscript indicates geodesic (no-SF) value]. For ro > 6M we have > 0, and the 
circular orbit is dynamically stable under small-e perturbations; the condition ujq = 
identifies the ISCO at ro — 6A/1^ To avoid confusion, it should be understood that 
the radial frequency too is a characteristic of the circular orbit (not the perturbed 
orbit), which, however, does not manifest itself in the dynamics of the circular orbit 
itself. Mathematically, ujq is associated with the curvature of the effective potential at 
its minimum; physically, it describes the radial frequency of an eccentric orbit as the 
eccentricity tends to zero. In our current discussion, we simply use the value of ujq as 
a convenient handle on the location of the ISCO. 

Now consider the effect of the SF, ignoring its dissipative piece. Equation p09p 
becomes 

^ir^T,s{r,Lz)+F:^^,. (112) 

The quantity Lz {— by its definition) is no longer constant along the orbit; its 
(conservative) time-variation is determined by the Lp component of the equation of 
motion, reading 

IiLz=F;°^\ (113) 

To learn how the SF affects the radial frequency of slightly eccentric orbits, we need 
to consider the linear variation of Eqs. (|112p and (|113p with respect to e. Through 
0(e), the relevant SF components along the orbit assume the general form 

Kons - ^o" + cosLOoT, = euoF^ sinwoT (114) 

§ The fact that Eq. gives uJq > also in the range ro < 3M is irrelevant here, since there are 

no timelike circular orbits in that range. 
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[this is shown in Ref. based on Eq. ([H])], where the coefhcients Fq, F[ and 
depend only on tq. In these expressions we can use the unperturbed radial frequency 
ojq instead of the perturbed frequency since the F coefficients are already 0{fi^), and 
it is only the leading-order SF effect which concerns us here. The linear variation 
procedure once again yields an equation of the form (|110p where now the radial 
frequency, perturbed by the conservative SF, is found to be 

Lu^ = ujI + {aF^ + f3F[ + ^F^) , (115) 

with a = -SrQ^ro - 4M)/(ro - 3M), /? = r^^ and 7 = -2ra^[M{ro - 3M)]i/2. This 
formula describes the 0(/i) conservative shift in the radial frequency off its geodesic 
value. Note it requires knowledge of the SF through 0(e). 

The perturbed ISCO radius, tq = 6A/+ Arjsco, is now obtained from the condition 
oj{ro) — 0. Recalling Eqs. ()111|) . this gives 

An,co = fi-\rl/M){3M ~ ro){aFS + f3F[ + lF^)l^^,^, 

= fi-^ (216M^FS -lOSM^F^ + V3fA , (116) 

V ^/ ro=6M 

where the substitution j'q — 6M is allowed since this only introduces an error of 
0{fi^) on the right-hand side. This equation describes the 0(/i) shift in the location 
of the ISCO due to the effect of the conservative SF. This shift is well defined 
albeit gauge dependent. However, the value of the azimuthal orbital frequency 
at the (shifted) location of the ISCO is gauge invariant. To obtain the shift in fl we 
simply substitute vq = 6M -f Arjsco in Eq. (|102p . writing = rio,isco + Arijsco where 
f2o,isco = na{6M) = (63/2 g^^^ through O(^), 

= -iAn_/M - ^M^,-^F,]i6M). (117) 

In Ref. the SF coefficient Fq , F[ and F^ are extracted numerically using the 
Lorenz-gauge SF code reviewed here in Sec. [Si With these numerical results, Eqs. 
(fTT6)) and (fTTT)) give the final values 

Arisco = -3.269 m, ^ 0A870 h/M. (118) 

^ 'O,isco 

We remind that the value of Arisco is specific to the Lorenz gauge, but that for Afiisco 
is invariant under all gauge transformations related to the Lorenz gauge through a 
transformation satisfying Eq. (|10ip . 

It is interesting to compare the conservative shift Afiisco with the frequency 
bandwidth of the dissipative transition across the ISCO, calculated by Ori and 
Thorne in [124] . Denoting the latter (in the Schwarzschild case) by Afidiss, one finds 
Afidiss/Afiisco ~ 9 (Ai/M)-3/5, giving, for example, - 35830, 9000 and 2261 for mass 
ratios n/M = 10^^, 10^^ and 10~^, respectively. This confirms our expectation that 
the conservative shift in the ISCO is practically negligible from the observational 
point of view. The main practical value of the results in Eq. (|118p remains that they 
provide an accurate strong-field benchmark to inform the development of approximate 
methods. 
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9. Reflections and prospects 

We have attempted here to capture a snapshot of the activity surrounding the 
development of rehable, efficient and accurate computational methods for the 
gravitational SF in black hole spacetimes. The problem still attracts considerable 
attention (more than half the items on our bibliographic list date > 2005), with a 
multitude of different approaches being studied by different groups. This multitude 
offers the opportunity for cross-validation of techniques and results — a particularly 
important prospect given the intricate nature of the numerics involved and the fact 
that SF calculations explore a new territory in black hole physics, yet uncharted 
neither by PN theory nor by Numerical Relativity. 

Indeed, the field has by now matured sufficiently that such cross-validation 
exercises are becoming possible. As we described in Sec. [H the last year has seen 
first quantitative comparisons between results from different calculations carried out 
in different gauges and using different numerical methods. We are now able to use SF 
codes to explore, for the first time, the conservative post-geodesic dynamics of strong- 
field orbits around a Schwarzschild black hole [53l [54l 11471 [26] . We can compute 
physical gauge-invariant SF effects and test them directly against results from PN 
theory in the weak- field regime [53l [Ml 11471 llOOj . Indeed, we can now start to use 
strong-field SF results in order to calibrate PN methods and assess their performance 
[25] . The exciting prospects for synergistic interaction between SF and PN theories are 
beginning to materialize, with much scope for important developments in the coming 
years. 

At preset, the state of the art in SF calculations is a code that can calculate the 
gravitational SF along any bound geodesic of the Schwarzschild geometry (currently 
at substantial computational cost, which future developments in the numerical 
technology may help reduce). This code, as many others mentioned in our review, 
is an implementation of the mode-sum regularization method (Sec. [4]), which has 
proven a useful framework for calculations in Schwarzschild. 

The Kerr problem, on the other hand, has hardly been tackled so far, and it 
represents the next significant challenge. Although the standard mode-sum method 
is in principle applicable to the gravitational SF in Kerr spacetime, the details of its 
numerical implementation in this case are yet to be worked out. It is possible that 
higher-dimensional techniques (like the m-mode scheme discussed in Sec. [7]) could 
provide an attractive alternative to standard mode-sum in the Kerr case. 

In the short term, activity is likely to focus on the following tasks: (i) continue 
to improve the computational efficiency of SF calculations using advanced numerical 
techniques (mesh refinement, spectral methods, etc.); (ii) tackle the Kerr problem; (iii) 
use SF codes to study post-geodesic physical effects (such as the finite-mass correction 
to the orbital precession rate), and in particular assess the relative importance of 
conservative SF effects in the EMRI problem; (iv) explore what can be learned from 
a comparison between SF and PN results. 

Within the wider context of the LISA EMRI problem, the computation of the 
SF on momentarily-geodesic particles is only one essential ingredient. There is still 
much more to understand before a faithful model of an astrophysical inspiral can 
be developed. Most crucially, a reliable and workable method must be devised for 
calculating the long-term evolution of the inspiral orbit. Work to address this problem 
has started recently [1361 11371 [571 [77] but much further development is required. 
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Appendix: Derivation of the regularization parameters 

We describe here the derivation of the regularization parameters for the gravitational 
SF in Kerr. The values of these parameters [stated in SeclH Eqs. ([Ml), dSS]) and gl])] 
were first published in Ref. [24j , but the detailed derivation has not appeared in print 
so far. We reproduce it here using the original method of Ref. [24]. 

Let the arbitrary timelike geodesic P be the worldline of a particle with mass n 
in a Kerr geometry with mass M ^ fj, and arbitrary spin aM. We wish to calculate 
the regularization parameters for the gravitational SF acting on the particle at an 
arbitrary point z along F with Boyer-Lindquist coordinates {tQ,ro,0o,(po). Without 
loss of generality we shall take (po = 0. 

We remind that the parameters ±^4", and C" are the leading-order coefficients 
in the formal expansion of the i-modes F^^J^^I^z) at large I [recall Eq. We also 

remind that the difference F^"j.^(z) — Fgj.{z) is expected to die off at large I faster 
than any power of l/l [recall the discussion surrounding Eq. ([55)) ]. Therefore, the 
values of A" , and C" can also — more conveniently — be deduced from the large-^ 
asymptotics of the S field modes: 

F^liz) = ±A°'L + 3" + C"/L + C»(L-2). (A.l) 

Once A", and are known, the parameter Z?" is given as the residual quantity 
[Eq. dST])] 

OC 

D" = ^[F^I^{z)tLA°' ~ B°' -C^]. (A.2) 

(=0 

Our starting point will be the local expression for the S-field h^p, given in Eq. 
(|19p. Referring back to Fig. [U we denote the Boyer-Lindquist coordinate difference 
between point x (an arbitrary point in the immediate vicinity of z) and point z by 
6x°' = x°' ~ z°'. In Eq. ([TO]) the quantities e^{x;z) and Ua{x;z) are smooth functions 
of (5x, and we may expand them in the form 

= So + Si+0{5x^), u„ = w„ + 0(fa2), (A.3) 

where 5*0 and 5*1 are, respectively, quadratic and cubic in 5x, and 5ua is linear in 5x. 
Explicitly, these expansion terms read 

5*0 ={gc.p + Uo,up)5x''5x^ = P^05x'^5xP, (A.4) 

= TlpPx^5x°'5xl^5x^ = Pc,p^Sx°'5x'^5x'^, (A.5) 

5u„ = r^pUxSx'^, (A.6) 
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where the background metric and connections are evaluated at z, and the coefhcients 
Pai3 and Pa/3j are those defined in Eq. (|^ (for a derivation of 5*1 see, for example, 
Appendix A of [55]). 

We now substitute the above expansions (|A.3p in Eq. (fTOj) . and subsequently 
construct the field Fg{x) = ^J''^x^^'h^{x), where, recall, the operator V'^^'^ is the 
one defined immediately below Eq. ([23]) . The resulting expression for Fg{x) can be 
written down as a sum of terms sorted according to how they scale with Sx: 

P^uiSx) Pu\{5x) PmiSx) 

^0 ^0 

Here eo = 'S'o j and P|-"-) denotes a certain multilinear function of the coordinate 
differences fe", of homogeneous order n in 6x. Note that the first term on the right- 
hand side scales as the second as 5x~^ , and the third as The 0{5x) 
remainder disappears at the limit x ^ z and cannot affect the value of the final SF; 
it is therefore safe to ignore it. The explicit form of P^"^ will not be needed in our 
analysis. The two other coefficients read 

P("i) ^ -2uf3U^Vf-'So = -Pyx^ (A.8) 

and 

(4) 



= \ [Pli'iPpPi + '^PppPi) " P'"\2Pxp[3 + Ppf}x)P^s] 5xP5xHx-^5x^. (A.9) 

The second equality in each of (jA.8[) and (|A.9[) was evaluated with the help of the 
following identities, which are easily confirmed: 

u^u^Vf = 7P"'V5, P^^Pp^ = P" , (A.IO) 



{5uiiu^ + Uf}5u^)V'^'^'^ Sq = -fc^u'^P'^Ja;'' --P'^pP-ySx'^Sx^ , (A.12) 



V^S-o = 2P5^feT, V55i - (2P5„/3 + P„;35)5a;"fe'', (A.ll) 

1 

2' 

where P^ is given in Eq. (|43p . 

To obtain the regularization parameters, we need to consider the decomposition 
of F^{x) in spherical harmonics (and then take x z). Notwithstanding the 
spheroidicity of the Kerr background, the spherical harmonics Yi,n{9, ip) are defined as 
usual on surfaces of constant r and t, where {t, r, 6, ip) are Boyer-Lindquist coordinates. 
To make this ^-mode decomposition easier, we introduce new coordinates: In general 
(i.e., for general Oq), the multipole decomposition will involve all azimuthal numbers 
\m\ < I for each I. However, if we choose a new set of angular coordinates, 
{0,(p) — > {9', If'), such that in the new coordinates z is positioned at the pole (i.e., 
0Q = 0), then the only contribution to each / mode in the new system would come from 
the axially-symmetric m = mode alone. Yet, due to the invariance of the Legendre 
functions under rotations on the 2-sphere, the overall / mode contribution (summed 
over m) would be the same in both systems. A transformation {0,ip) — > {9', if') is 
given explicitly by 

cos 9' — cos ip sin 9 sin + cos 9 cos 6*0 , 
, , sin ip sin 9 
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where ip'^ indicates the azimuthal direction of the particle's velocity in the new polar 
system (i.e., in the new system the velocity is momentarily tangent to the longitude 
line (fi' = const = (p'^ on the 2-sphere). The angle tp'^ depends on and u"^, but the 
explicit relation will not be needed here. 

Since the new angular coordinates are not well behaved at z, we also introduce 
Cartesian-like local coordinates based around z: 

X = p{e') cos((p' - (^;), y = p{e') sm{p' - ^'J, (A.14) 

where we require p — 9' + 0{9'^) near z. Note that at z we have a; = y = as well as 
— 0. For later convenience, we make here the concrete choice 

p^2sm{6'/2), (A.15) 

giving 

cos6i' = 1 - i(a;2 +2/2). (A.16) 

In the following we will need the transformation [9, ip) ^ {x, y) only in the 
immediate neighborhood of z. To sufficient order, we find 

69 =x+ (1/2) cot6'oy2 + 0{5x^), 

Sip = (sin 9oy\y - cot Boxy) + 0{5x^), (A.f7) 

where, recall, 59 = 9 — 9q and 5'p — Lp — (po = ip are the Boyer-Lindquist coordinate 
differences, and 0{5x^) represents terms at least cubic in x and y. We denote 

X" = (X",X\X2,X3) = {5t,5r,x,y/ smOo), (A.18) 

noting X" = at the particle. The direct transformation from X" to the Boyer- 
Lindquist coordinate differences (5a;" can be written to sufhcient order in the compact 
form 

5x°' = X" + C;^^X^X^ (A. 19) 

where the constant coefficients are those defined in Eq. (|^^ . 

Our next step is to express the S-field force Fc^{x) in terms of the X°' coordinates 
(note, however, that we do not consider here the X" components of Fg'(x), but rather 
the Boyer-Lindquist components; the X°' coordinates are merely introduced to assist 
with the Legendre integration below) . We substitute Eq. ()A.I9p in Eq. (IA.7P and then 
re-expand in the coordinate differences X" . The result takes the form 



where 



ii-^Fi{x) =FX + FS+FS+ OiSx), (A.20) 



pa x^' 

n - — (A.2f) 



P''„a...XPX'^X^X^ 



F§ = ^ P^^^^^;^ , (A.22) 

P",n,x.R^f,Xi'X''X^X''Xl^X^X^ 
Fg = ^"'^"'^^^ 7 . (A.23) 

^QX 



Here 



= Po^pX^X^, (A.24) 
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P^^^Xa-p-yS constant coefHcients whose explicit values will not be needed, and 

+ {iP^Ppx ~ P\Ppp)C^s- (A.25) 
[In the last expression, the first line is simply the coefficient of P"^-^ from Eq. (jA.9[) . 
and the second line arises from expanding the term t^^P^^^^ in Eq. (|A.7p in powers of 
X".] Note in Eq. (|X20)) that the terms F%, and Fg scale as X~'^, and 
X~^, respectively. The terms included in 0{Sx) vanish at x ^ z, and are therefore 
irrelevant for calculating the regularization parameters — we can hence safely ignore 
them in what follows. 

To obtain the regularization parameters, via Eqs. (jA.ip and (|A.2p . we now need to 
consider the spherical-harmonic I modes F^j_{z). These are constructed, by definition, 
through 

Fgl= lim — / d cos 61' / dip'Pi{cos9')F^{r,9',(p';z), (A.26) 

where P;(cos6'') is the Legendre polynomial and, recall, L = I + 1/2. Here we have 
already taken the limits t — > tg, 9 ~^ 9q and (/? — > (^q? thereby choosing to approach 
the particle along the radial direction, recalling that we expect two different one-side 
values (hence the label '±'). Let us write 

^-^Fil = Ff^ + < + FS\ (A.27) 

where the three terms on the right-hand side represent the respective contributions 
to F^l_ from the three terms F^ ^ q in Eq. ((X20I) . As explained in Ref. the 
second and third terms (scaling near the particle as ^ Sx^^ and ^ 6x'^ , respectively) 
are sufficiently regular to allow interchanging the order of the limit and integration 
m Eq. 1X2611 . For the same reason, the radial limit is well defined and the two-side 
ambiguity does not occur. As we discuss below, such ambiguity does appear when 
considering F^^^ — hence the label ±. In what follows we consider each of the three 
contributions to F^j_ in turn. 

We start with Fg'^ . This term is obtained by replacing F^ in Eq. (|A.26p with Fg 
from Eq. (jA.23p . and we are allowed to pull the limit 6r 0^ over the integrals. The 
outcome has the form 

^c^^l ' dcos9' l\^'Pi{co8 9')e~l{x,y)P('^){x,y), (A.28) 

where P^^{x, y) is a certain polynomial of homogeneous order 7 in the two "angular" 
coordinates {x,y) = (A"^, AT'^ sin6'o), and 

eox = eox(<5r ^ 0, (5i -> 0) = [[g^^{z) + ul\x^ ^ 9yyiz)y^] ^'^ ■ (A.29) 

In the last equality we used Eq. (|A.24[) with Pa^ — gap + UaUp, recalling Uj^ = 
along with gxy{z) = 0. We observe that egx is an even function of both x and y, 
and, recalling Eq. (|A.16p . so is cos6''. However, each of the terms in P^iAx,y) (such 
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as oc x'^y^, cx x'^y^, etc.) is odd in either x or y. It then immediately foUows from 
symmetry that the integral in Eq. (|A.28p vanishes: 



Fg^ = 0. 



(A.30) 

and (5t — > 0, we have from Eq. 



Next consider F^^. First, taking the limits 6r 
(1X221) 

F^iSr, St^O)^ i-^P^.^.X^X^X^X", (A.31) 

where roman indices run over the angular coordinates {X'^,X^) only. Then, using 
{X^,X^) — (x, y/ sin 6*0) = /ci(cos 93', siniys'/ sin^o) = pw"", we write X" = pw°' and 
^ox ~ P^PabW°"w^ ■ This allows us to write the Legendre integral for Fg^ in a factorized 
form: 



— dcos0 7— — 

27r A p{e') 



dtp 



/ ^ a 



bed 



{PabW^W^f/^ 



(A.32) 



Here the 0' integral is elementary: The entire term in the first set of square brackets 
reads simply {2tt)~^ (for any /). The term in the second set of square brackets is 



PZh^J"'""^, where I"'""^ are the (^-independent) integrals given in Eq. These 
integrals, recall, are not elementary, but they can be expressed explicitly in terms of 
complete elliptic integrals, as in Eq. (|49| . In conclusion, we find 



pal _ (n \-lpa jw 
— [ATT) r- abcd^ 



bed 



(A.33) 



Importantly, the term F^^ contributes to each Z-mode F^^ a constant amount, 
independent of I. 

We have one more contribution to Fg*' to consider: that of -F^i- Recalling Eq. 
(|A.2ip . we have 

^A± = j™^ ^ / Pi (^os ^')eoin^^ (A-34) 



where the integral is taken over the 2-sphere, we have used the Jacobian 
d{cos 9' ,ip')/d{x,y) — —1, and the integrand is understood to be already evaluated at 
t = to. The Legendre polynomial can be written in the form P;(cos0') = 1 + p^H{p), 
where H{p) admits a regular Taylor expansion in p^. Consider first the contribution 
of the term c>c p^H{p) to -Fa±' corresponding integrand in Eq. (|A.34p has the 
form e'^^H{p)P"^A5r,x,y) (recalling p^ = x"^ + y^), where P^^^ is a polynomial of 
homogeneous order 3 in 5r, x, y. This contribution to the integrand is therefore 
bounded, and we are allowed to swap the limit and integral just as we did with 
FqI^ . The resulting integral then vanishes by virtue of the same symmetry argument 
applied in the case of -Fq': Both functions e^^ and H{p) are even in each of x and y, 
while all the possible terms in Pl^^{0, x,y) are odd in either x or y. 
We are thus left with the contribution 



Fa±= lini 



L 



dxdyeQ^P°:.X^'. 



(A.35) 



To evaluate this integral, we divide the integration domain (the 2-sphere) into two 
regions: Let Dm denote the square —h < x,y < h for some particular < /i < 1 
(say, h = 1/10); and let Uout denote the remaining integration area. Consider first 
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the contribution to from -Dout ■ Since in this domain the integrand is regular (the 
only singularity is at x = y = 0, which is in I?in), we are allowed to interchange the 
limit and integration. As a result, the integrand takes the form e^^P^X'', and the 
integral over -Dout vanishes by virtue of the odd symmetry. The remaining piece of 
the integral, over Z?in, is 

^A± - , lini ^ f f dx dy e-^P^X^^. (A.36) 

This integral takes a simpler form if we express it in terms of the rescaled coordinates 
X = x/Sr, y = y/Sr, and A"^ = X^/Sr = {0,l,x,y/ s'mOo) (where we have already 
taken St 0). For given (5r 7^ we have eox = {PapX^'X'^y/^ = ±Sr{Paf3X"'X'^)^/^ , 
where the sign corresponds to the sign of Sr. Hence we obtain 

Fa± = ± , lim^ — / / didy ^ - 



5r^0± 2tT j_,^/sr J-h/Sr {PapX^XI')^^ 

L r°° p°'„xi^ 

±— / didy - , , (A.37) 

where the last equality holds because the integrand, expressed in terms of the tilde 
variables, no longer depends on Sr. This integral is now elementary, and evaluating it 
gives 

F^^t = ±LA°', (A.38) 

where the signs correspond to Sr — > 0^, and where the various Boyer-Lindquist 
components of A" are given in Eq. ([55)1 . [Section V.D of Ref. [55] explains in detail 
how the integral in Eq. (jA.37[) is evaluated in the special case of Schwarzschild, and 
the method is directly applicable to Kerr.] Notice that is found to depend on I 
only through the prefactor L. 

In summary, collecting the results ([Ool . ([Al30| . (|Al33| and (jXaS]) . we find that 
the I modes Fg'^ are given precisely by 

F^l = ±LA°' + B", (A.39) 

where A" and B" = {2Tr)-^ P^^^I"'""^ are /-independent. Comparing with Eq. (|A?T|) . 
we identify A" and as the first two of the regularization parameters. This 
comparison also implies C" = 0. Furthermore, we find that each of the individual 
terms in the sum over I in Eq. (|A.2p vanishes, giving also = 0. This completes the 
derivation of all regularization parameters. 

We comment on a potentially confusing aspect of the above analysis: We have 
discarded in Eq. (jA.20p the 0{Sx) terms of Fg'_^_{x), which vanish at x ^ z. Clearly, 
the multipole expansion of these neglected terms could contribute to F^j. [e.g., they 
may well add a term cx L^^ in Eq. (jA.39|l ]. Such terms, however, must add up to zero 
upon summation over I (when evaluated at the particle). They hence affect neither 
the value of Z?" in Eq. (K^ . nor the value of the final SF in Eq. (pSD . 
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